Commit 42d3d8a2 authored by Robin Brost's avatar Robin Brost
Browse files

working on the evaluation method

parent 9f152982
Pipeline #12499 failed with stages
in 7 minutes and 41 seconds
<?xml version="1.0" encoding="UTF-8"?>
<!DOCTYPE QtCreatorProject>
<!-- Written by QtCreator 4.9.1, 2019-10-17T17:21:30. -->
<!-- Written by QtCreator 4.9.1, 2019-10-21T16:15:44. -->
<qtcreator>
<data>
<variable>EnvironmentId</variable>
......
......@@ -106,15 +106,22 @@ int main(void)
b.setZero();
// first constraint: first variable equals three times second
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;
//different number of constraints :
if(n_constraints == 1){
A.coeffRef(0,0) = 2;
A.coeffRef(0,1) = -6;
b.coeffRef(0) = 0;
}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(2,1) = -3;
b.coeffRef(0) = 0;
b.coeffRef(1) = 27;
b.coeffRef(2) = 0;
}
std::cout << A << std::endl << b << std::endl;
......@@ -144,10 +151,13 @@ int main(void)
Eigen::VectorXi d = b.cast<int>();
ExactConstraintSatisfaction satisfy;
satisfy.printMatrix(C);
std::cout << std::endl;
satisfy.IREF_Gaussian(&C,&d);
satisfy.printMatrix(C);
std::cout << std::endl;
satisfy.IRREF_Jordan(&C,&d);
satisfy.printMatrix(C);
std::cout << std::endl;
satisfy.printVector(d);
x.coeffRef(0) = 3.0 * x.coeffRef(1);
......
......@@ -57,6 +57,49 @@ void ExactConstraintSatisfaction::printVector(Eigen::VectorXi b)
std::cout << std::endl;
}
int ExactConstraintSatisfaction::largestExponent(Eigen::VectorXd x)
{
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);
}
largest_exponent = expo;
delta = std::pow(2, expo);
return expo;
}
double ExactConstraintSatisfaction::F_delta(double x)
{
int sign = -1;
if(x >= 0)
sign = 1;
double x_of_F = (x + sign * delta) - sign * delta;
return x_of_F;
}
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 lcm_D = 1;
for(int d : D){
lcm_D = lcm(lcm_D, d);
}
return lcm_D;
}
int ExactConstraintSatisfaction::indexPivot(Eigen::SparseMatrix<int>::RowXpr row)
{
for(int i = 0; i < row.size(); i++){
if(row.coeffRef(i) != 0)
return i;
}
return -1;
}
// ----------------------Matrix transformation-----------------------------------
void ExactConstraintSatisfaction::IREF_Gaussian(Eigen::SparseMatrix<int>* A, Eigen::VectorXi* b){
number_pivots = 0;
......@@ -105,7 +148,6 @@ void ExactConstraintSatisfaction::IREF_Gaussian(Eigen::SparseMatrix<int>* A, Eig
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;
......@@ -132,93 +174,78 @@ void ExactConstraintSatisfaction::IRREF_Jordan(Eigen::SparseMatrix<int> *A, Eige
}
}
void ExactConstraintSatisfaction::evaluation(Eigen::SparseMatrix<int> *A, Eigen::VectorXi *b, Eigen::VectorXi x)
{
//-------------------------------------Evaluation--------------------------------------------------------------------
}
void ExactConstraintSatisfaction::makeDiv(std::set<int> D, double x)
void ExactConstraintSatisfaction::evaluation(Eigen::SparseMatrix<int> *A, Eigen::VectorXi *b, Eigen::VectorXd x)
{
int cols = A->cols();
largestExponent(x);
for(int k = cols -1; k >= 0; k--){
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){
pivot = i; //the row with the pivot
break;
}
void ExactConstraintSatisfaction::safeDot(std::set<ExactConstraintSatisfaction::pair> S)
{
}
if(pivot == -1){ //there is no pivot in this column
std::set<int> D;
D.clear();
for(int i = 0; i <= std::min(k, number_pivots - 1); i ++){
if(A->coeffRef(i,k) != 0){
D.insert(A->coeffRef(i, indexPivot(A->row(i))));
}
}
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;
S.clear();
for(int i = k+1; i < cols; i++){ //construct the Set S to do the dot Product
std::pair<double, double> tuple;
tuple.first = A->coeffRef(pivot,i);
tuple.second = x.coeffRef(i) / A->coeffRef(pivot,k);
S.insert(tuple);
}
x.coeffRef(k) = b->coeffRef(k) - safeDot(S);
}
}
}
double ExactConstraintSatisfaction::makeDiv(const std::set<int>& D, double x)
{
if(D.empty()){
return F_delta(x);
}
int d = lcm_Set(D);
return F_delta(x/d) * d;
/*
// ------------------------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;
double ExactConstraintSatisfaction::safeDot(const std::set<std::pair<double, double>> &S)
{
std::set<std::pair<double, double>> P, N;
P.clear();
N.clear();
for(std::pair<double, double> element : S){
if(element.first * element.second > 0)
{
P.insert(element);
}
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));
else
{
N.insert(element);
}
gcdValue = gcd(gcdValue, b);
return gcdValue;
}
}
double r = 0;
while(!P.empty() || !N.empty()){
if(!P.empty() && (r < 0 || N.empty())){
//std::pair<double, double> element = P.begin();
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;
}
}
}
*/
......@@ -2,12 +2,13 @@
#define EXACTCONSTRAINTSATISFACTION_HH
#include <CoMISo/Config/config.hh>
#include <CoMISo/Config/CoMISoDefines.hh>
#include <CoMISo/NSolver/NProblemInterface.hh>
#include <vector>
#include <set>
class ExactConstraintSatisfaction
class COMISODLLEXPORT ExactConstraintSatisfaction
{
public:
ExactConstraintSatisfaction();
......@@ -15,10 +16,8 @@ public:
//-----------------------helpfull variables-------------------------------
int number_pivots = 0; //number of rows with a pivot;
struct pair{
int x;
int y;
};
int largest_exponent = 0;
double delta = 0;
//-----------------------helpfull methods---------------------------------
......@@ -31,6 +30,12 @@ public:
void printMatrix(Eigen::SparseMatrix<int> A);
void printVector(Eigen::VectorXi b);
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 indexPivot(Eigen::SparseMatrix<int>::RowXpr row);
//--------------------matrix transformation-------------------------------
void IREF_Gaussian(Eigen::SparseMatrix<int>* A, Eigen::VectorXi* b);
......@@ -38,9 +43,9 @@ public:
//-------------------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);
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);
};
#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