preliminaries for variable restricted dualization: Namely Coefficients with fixed Variables

Linear_Programming/Coefficient.cpp

+#include "Coefficient.h"
+
+Coefficient::Coefficient(std::map<Monomial, double, Monomial::less> coeff) : _coeff(coeff) {}
+
+Coefficient::Coefficient(){}
+Coefficient::Coefficient(double initial_val) : _coeff({{{}, initial_val}}) {}
+
+
+bool Coefficient::is_non_zero(){
+  for(std::pair<const Monomial, double>& m : this->_coeff){
+    if(std::abs(m.second) > 1e-100){
+      return true;
+    }
+  }
+  return false;
+}
+
+bool Coefficient::is_unity(){
+  auto search_result = this->_coeff.find({});
+  if(search_result != this->_coeff.end()){
+    if(std::abs(search_result->second - 1) > 1e-100){
+      return false;
+    }else{
+      auto tmp = search_result->second;
+      search_result->second = 0;
+      bool non_unity = this->is_non_zero();
+      search_result->second = tmp;
+      if(non_unity){
+        return false;
+      }
+    }
+  }
+  return true;
+}
+
+// Coefficient operator+(double val, Coefficient coeff){
+//   Coefficient tmp (coeff);
+//
+//   for(std::pair<Monomial, double>& m : tmp){
+//     m.second += val;
+//   }
+//
+//   return tmp;
+// }
+Coefficient operator+(Coefficient coeffa, Coefficient coeffb){
+  Coefficient tmp (coeffa);
+
+  for(std::pair<const Monomial, double>& m : coeffb._coeff){
+    auto search_result = tmp._coeff.find(m.first);
+    if(search_result != coeffb._coeff.end()){
+      search_result->second += m.second;
+    }else{
+      tmp._coeff.insert(m);
+    }
+  }
+
+  return tmp;
+}
+Coefficient operator*(double val, Coefficient coeff){
+  Coefficient tmp (coeff);
+
+  for(std::pair<const Monomial, double>& m : tmp._coeff){
+    m.second *= val;
+  }
+  return tmp;
+}
+Coefficient operator*(Coefficient coeffa, Coefficient coeffb){
+  Coefficient tmp (coeffa);
+
+  for(std::pair<const Monomial, double>& m : coeffb._coeff){
+    for(std::pair<const Monomial, double>& n : coeffa._coeff){
+      tmp._coeff.insert({m.first*n.first, m.second*n.second});
+    }
+  }
+
+  return tmp;
+}
+
+//Coefficient operator-(double val, Coefficient coeff);
+Coefficient operator-(Coefficient coeffa, Coefficient coeffb){
+  Coefficient tmp (coeffa);
+
+  for(std::pair<const Monomial, double>& m : coeffb._coeff){
+    auto search_result = tmp._coeff.find(m.first);
+    if(search_result != coeffb._coeff.end()){
+      search_result->second -= m.second;
+    }else{
+      tmp._coeff.insert(m);
+    }
+  }
+
+  return tmp;
+}
+
+Coefficient operator-(Coefficient coeff){
+  Coefficient tmp;
+  for(std::pair<const Monomial, double> m : coeff._coeff){
+    m.second = -m.second;
+    tmp._coeff.insert(m);
+  }
+
+  return tmp;
+}
+
+std::ostream& operator<<(std::ostream& os, Coefficient coeff){
+  if(coeff.is_non_zero()){
+    bool parenthesis = coeff._coeff.size() > 1 /*|| (coeff._coeff.size() > 0 && coeff._coeff.begin()->second < 1e-100)*/;
+    if(parenthesis){
+      os << "(";
+    }
+    bool add_plus = Coefficient::always_add_sign;
+    for(std::pair<const Monomial, double>& m : coeff._coeff){
+      if(std::abs(m.second) > 1e-100){
+        if(m.second < 1e-100){
+          if(add_plus){
+            os << " - ";
+          }else{
+            os << "-";
+          }
+        }else if(add_plus){
+          os << " + ";
+        }
+        add_plus = true;
+
+        if(std::abs(m.second-1) < 1e-100){
+          os << m.first;
+        }else{
+          os << std::abs(m.second) << " " << m.first;
+        }
+      }
+    }
+    if(parenthesis){
+      os << ")";
+    }
+  }
+
+  return os;
+}

Linear_Programming/Coefficient.h

+#ifndef COEFFICIENT_H
+#define COEFFICIENT_H
+
+#include <cstddef>
+#include <map>
+#include <iostream>
+#include <string>
+#include <cassert>
+#include <stdexcept>
+#include <cmath>
+
+#include "Monomial.h"
+#include "Variable.h"
+
+class Coefficient{
+  std::map<Monomial, double, Monomial::less> _coeff;
+
+  Coefficient(std::map<Monomial, double, Monomial::less> coeff);
+public:
+  friend Coefficient operator+(double val, Coefficient coeff);
+  friend Coefficient operator+(Coefficient coeffa, Coefficient coeffb);
+  friend Coefficient operator*(double val, Coefficient coef);
+  friend Coefficient operator*(Coefficient ceoffa, Coefficient coeffb);
+
+  friend Coefficient operator-(double val, Coefficient coeff);
+  friend Coefficient operator-(Coefficient coeffa, Coefficient coeffb);
+  friend Coefficient operator-(Coefficient);
+
+  friend std::ostream& operator<<(std::ostream& os, Coefficient coeff);
+
+  Coefficient();
+  Coefficient(double initial_val);
+
+  bool is_non_zero();
+  bool is_unity();
+  double value();
+
+  static const bool always_add_sign = true;
+};
+
+std::ostream& operator<<(std::ostream& os, Coefficient coeff);
+
+#endif

Linear_Programming/Constraint.cpp

 #include "Constraint.h"
 
-Constraint::Constraint(std::string description, relation relation_type, std::map<Variable*, double> lhs, double rhs)
+Constraint::Constraint(std::string description, relation relation_type, std::map<Variable*, Coefficient> lhs, Coefficient rhs)
   : _description(description), _relation(relation_type), _lhs(lhs), _rhs(rhs){}
 
-Constraint::Constraint(relation relation_type, std::map<Variable*, double> lhs, double rhs)
+Constraint::Constraint(relation relation_type, std::map<Variable*, Coefficient> lhs, Coefficient rhs)
   : _relation(relation_type), _lhs(lhs), _rhs(rhs){}
 
 std::string Constraint::description(){
@@ -14,19 +14,19 @@ relation Constraint::is_equality(){
   return this->_relation;
 }
 
-std::map<Variable*, double> Constraint::lhs(){
+std::map<Variable*, Coefficient> Constraint::lhs(){
   return this->_lhs;
 }
 
-double Constraint::lhs_coefficient(Variable* var){
+Coefficient Constraint::lhs_coefficient(Variable* var){
   auto tmp = this->_lhs.find(var);
   if(tmp == this->_lhs.end()){
-    return 0;
+    return {0};
   }
   return tmp->second;
 }
 
-double Constraint::rhs(){
+Coefficient Constraint::rhs(){
   return this->_rhs;
 }
 
@@ -36,18 +36,13 @@ std::ostream& operator<<(std::ostream& os, Constraint& constraint){
 
   bool add_plus = false;
   for(auto summand : constraint.lhs()){
-    if(std::abs(summand.second) > 1e-100){
+    if(summand.second.is_non_zero()){
       empty = false;
-      if(summand.second < 0){
-        os << " - ";
-      }else if(add_plus){
+      if(add_plus && !Coefficient::always_add_sign){
         os << " + ";
       }
       add_plus = true;
-      if(std::abs(std::abs(summand.second) - 1) > 1e-100){
-        os << std::abs(summand.second) << " ";
-      }
-      os << summand.first->description();
+      os << summand.second << *summand.first;
     }
   }
   if(empty){
@@ -58,6 +53,12 @@ std::ostream& operator<<(std::ostream& os, Constraint& constraint){
   }else{
     os << " <= ";
   }
-  os << constraint.rhs();
+  if(!constraint.rhs().is_non_zero()){
+    os << "0";
+  }else if(constraint.rhs().is_unity()){
+    os << "1";
+  }else{
+    os << constraint.rhs();
+  }
   return os;
 }

Linear_Programming/Constraint.h

@@ -8,23 +8,24 @@
 #include <cmath>
 
 #include "Variable.h"
+#include "Coefficient.h"
 
 enum relation : bool {Inequality = false, Equality = true};
 
 class Constraint{
   std::string _description;
   relation _relation;
-  std::map<Variable*, double> _lhs;
-  double _rhs;
+  std::map<Variable*, Coefficient> _lhs;
+  Coefficient _rhs;
 public:
-  Constraint(std::string description, relation relation_type, std::map<Variable*, double> lhs, double rhs);
-  Constraint(relation relation_type, std::map<Variable*, double> lhs, double rhs);
+  Constraint(std::string description, relation relation_type, std::map<Variable*, Coefficient> lhs, Coefficient rhs);
+  Constraint(relation relation_type, std::map<Variable*, Coefficient> lhs, Coefficient rhs);
 
   std::string description();
   relation is_equality();
-  std::map<Variable*, double> lhs();
-  double lhs_coefficient(Variable* var);
-  double rhs();
+  std::map<Variable*, Coefficient> lhs();
+  Coefficient lhs_coefficient(Variable* var);
+  Coefficient rhs();
 };
 
 std::ostream& operator<<(std::ostream& os, Constraint& constraint);

Linear_Programming/Linear_Program.cpp

@@ -10,24 +10,24 @@ Linear_Program::Linear_Program(Linear_Program&& lp){
 
 Linear_Program::Linear_Program(bool maximum) : _maximum(maximum) {}
 
-Linear_Program::Linear_Program(bool maximum, std::map<Variable*, double> direction, Polyeder polyeder)
+Linear_Program::Linear_Program(bool maximum, std::map<Variable*, Coefficient> direction, Polyeder polyeder)
   : _maximum(std::move(maximum)), _direction(std::move(direction)), _polyeder(std::move(polyeder)) {}
 
 bool Linear_Program::is_maximum(){
   return this->_maximum;
 }
 
-std::map<Variable*, double> Linear_Program::direction(){
+std::map<Variable*, Coefficient> Linear_Program::direction(){
   return this->_direction;
 }
 
-void Linear_Program::add_direction_coefficient(std::pair<Variable*, double> summand){
+void Linear_Program::add_direction_coefficient(std::pair<Variable*, Coefficient> summand){
   this->_direction.insert(summand);
 }
 
-double Linear_Program::direction_coefficient(Variable* var){
+Coefficient Linear_Program::direction_coefficient(Variable* var){
   auto search_result = this->_direction.find(var);
-  if(search_result == this->_direction.end()) return 0;
+  if(search_result == this->_direction.end()) return {0};
   return search_result->second;
 }
 
@@ -76,7 +76,7 @@ Linear_Program Linear_Program::relaxation_dual(){
     }
   } variables;
 
-  std::map<Variable*, double> direction;
+  std::map<Variable*, Coefficient> direction;
 
   for(size_t index = 0; index < this->polyeder().number_of_constraints(); ++index){
     std::stringstream name;
@@ -93,9 +93,9 @@ Linear_Program Linear_Program::relaxation_dual(){
   // we introduce one constraint per primal variable
   for(auto curr_variable : this->polyeder().variables()){
     // gather lhs_coefficient
-    std::map<Variable*, double> lhs;
+    std::map<Variable*, Coefficient> lhs;
     for(size_t index = 0; index < this->polyeder().number_of_constraints(); ++index){
-      lhs.insert({variables.unwrap().right.find(index)->second, -this->polyeder().constraint(index).lhs_coefficient(curr_variable.left)});
+      lhs.insert({variables.unwrap().right.find(index)->second, {-this->polyeder().constraint(index).lhs_coefficient(curr_variable.left)}});
     }
 
 
@@ -104,7 +104,6 @@ Linear_Program Linear_Program::relaxation_dual(){
     name << "dual_" << curr_variable.left->description();
 
     // generate contraint where inequality / equality property is based in existance of 0 <= x_curr_var bound
-
     {
       // if upper bound ex. generate new var and add to lhs
       if(curr_variable.left->upper_bound().first == true){
@@ -113,7 +112,7 @@ Linear_Program Linear_Program::relaxation_dual(){
         Variable* upper_bound_var = new Variable(ub_name.str(), Continuous, {true, 0});
         variables.insert(upper_bound_var);
 
-        lhs.insert({upper_bound_var, -1});
+        lhs.insert({upper_bound_var, {-1}});
         direction.insert({upper_bound_var, curr_variable.left->upper_bound().second});
       }
     }
@@ -128,9 +127,9 @@ Linear_Program Linear_Program::relaxation_dual(){
           Variable* lower_bound_var = new Variable(lb_name.str(), Continuous, {true, 0});
           variables.insert(lower_bound_var);
 
-          lhs.insert({lower_bound_var, 1});
+          lhs.insert({lower_bound_var, {1}});
           constraints.push_back(Constraint(name.str(), Equality, lhs, -this->direction_coefficient(curr_variable.left)));
-          direction.insert({lower_bound_var, -1});
+          direction.insert({lower_bound_var, {-1}});
         }
       }else{
         constraints.push_back(Constraint(name.str(), Equality, lhs, -this->direction_coefficient(curr_variable.left)));
@@ -157,19 +156,14 @@ std::ostream& operator<<(std::ostream& os, Linear_Program& linear_program){
   bool empty = true;
   if(linear_program.direction().size() > 0){
     bool add_plus = false;
-    for(auto direction_summand : linear_program.direction()){
-      if(std::abs(direction_summand.second) > 1e-100){
+    for(std::pair<Variable*, Coefficient> direction_summand : linear_program.direction()){
+      if(direction_summand.second.is_non_zero()){
         empty = false;
-        if(direction_summand.second < 0){
-          os << " - ";
-        }else if(add_plus){
+        if(add_plus && !Coefficient::always_add_sign){
           os << " + ";
         }
         add_plus = true;
-        if(std::abs(std::abs(direction_summand.second) - 1) > 1e-100){
-          os << std::abs(direction_summand.second) << " ";
-        }
-        os << direction_summand.first->description();
+        os << direction_summand.second << direction_summand.first->description();
       }
     }
   }

Linear_Programming/Linear_Program.h

@@ -17,18 +17,18 @@
 
 class Linear_Program{
   bool _maximum;
-  std::map<Variable*, double> _direction;
+  std::map<Variable*, Coefficient> _direction;
   Polyeder _polyeder;
 public:
   Linear_Program(Linear_Program& lp);
   Linear_Program(Linear_Program&& lp);
   Linear_Program(bool maximum);
-  Linear_Program(bool maximum, std::map<Variable*, double> direction, Polyeder polyeder);
+  Linear_Program(bool maximum, std::map<Variable*, Coefficient> direction, Polyeder polyeder);
 
   bool is_maximum();
-  std::map<Variable*, double> direction();
-  void add_direction_coefficient(std::pair<Variable*, double> summand);
-  double direction_coefficient(Variable* index);
+  std::map<Variable*, Coefficient> direction();
+  void add_direction_coefficient(std::pair<Variable*, Coefficient> summand);
+  Coefficient direction_coefficient(Variable* index);
   Polyeder& polyeder();
 
   void solve();
@@ -37,6 +37,8 @@ public:
   void operator=(Linear_Program&& lp);
 
   Linear_Program relaxation_dual();
+  Linear_Program Danzig_Wolfe_IPM_complete_optimal_inequality_relaxation();
+  Linear_Program Benders_Reformulierung_integral_relaxation();
 };
 
 std::ostream& operator<<(std::ostream& os, Linear_Program& linear_program);

Linear_Programming/Monomial.cpp

+#include "Monomial.h"
+
+Monomial::Monomial(const Monomial& other) : _vars(other._vars) {}
+
+Monomial::Monomial(std::set<Variable*> vars) : _vars(vars) {}
+
+bool Monomial::less::operator()(const Monomial a, const Monomial b) const {
+  std::stringstream key_a;
+  for(auto var_ptr : a._vars){
+    key_a << var_ptr;
+  }
+  std::stringstream key_b;
+  for(auto var_ptr : b._vars){
+    key_b << var_ptr;
+  }
+  return key_a.str() < key_b.str();
+}
+
+Monomial operator*(const Monomial a, const Monomial b) {
+  Monomial tmp (a);
+  tmp._vars.insert(b._vars.begin(), b._vars.end());
+  return tmp;
+}
+
+std::ostream& operator<<(std::ostream& os, const Monomial m){
+  for(Variable* v : m._vars){
+    os << v->description() << " ";
+  }
+
+  return os;
+}

Linear_Programming/Monomial.h

+#ifndef MONOMIAL_H
+#define MONOMIAL_H
+
+#include <set>
+#include <sstream>
+
+#include "Variable.h"
+
+class Monomial{
+  std::set<Variable*> _vars;
+public:
+  Monomial(const Monomial& other);
+  Monomial(std::set<Variable*> vars = {});
+
+  class less{
+  public:
+    bool operator()(const Monomial a, const Monomial b) const;
+  };
+  friend less;
+
+  friend Monomial operator*(const Monomial a, const Monomial b);
+  friend std::ostream& operator<<(std::ostream& os, const Monomial m);
+};
+
+
+#endif

Linear_Programming/Variable.cpp

@@ -3,13 +3,13 @@
 bool Variable::VERBOSE_IDENT = true;
 
 
-Variable::Variable(const Variable& var) : _description(var._description), _integrality(var._integrality), _lower_bound(var._lower_bound), _upper_bound(var._upper_bound), _value(var._value) {}
+Variable::Variable(const Variable& var) : _description(var._description), _integrality(var._integrality), _lower_bound(var._lower_bound), _upper_bound(var._upper_bound), _fixed(false), _value(var._value) {}
 
-Variable::Variable(Variable&& var) : _description(std::move(var._description)), _integrality(std::move(var._integrality)), _lower_bound(std::move(var._lower_bound)), _upper_bound(std::move(var._upper_bound)), _value(std::move(var._value)) {}
+Variable::Variable(Variable&& var) : _description(std::move(var._description)), _integrality(std::move(var._integrality)), _lower_bound(std::move(var._lower_bound)), _upper_bound(std::move(var._upper_bound)), _fixed(false), _value(std::move(var._value)) {}
 
 //Variable::Variable(integrality integrality, std::pair<bool, double> lower_bound, std::pair<bool, double> upper_bound, double value) : _integrality(integrality), _lower_bound(lower_bound), _upper_bound(upper_bound), _value(value) {}
 
-Variable::Variable(std::string description, integrality integrality, std::pair<bool, double> lower_bound, std::pair<bool, double> upper_bound, double value) : _description(description), _integrality(integrality), _lower_bound(lower_bound), _upper_bound(upper_bound), _value(value) {}
+Variable::Variable(std::string description, integrality integrality, std::pair<bool, double> lower_bound, std::pair<bool, double> upper_bound, double value) : _description(description), _integrality(integrality), _lower_bound(lower_bound), _upper_bound(upper_bound), _fixed(false), _value(value) {}
 
 
 std::string Variable::description(){
@@ -29,6 +29,20 @@ std::pair<bool, double> Variable::upper_bound(){
 }
 
 
+bool Variable::is_fixed(){
+  return this->_fixed;
+}
+
 double& Variable::value(){
+  assert(this->is_fixed());
   return this->_value;
 }
+
+std::ostream& operator<<(std::ostream& os, Variable& var){
+  if(var.is_fixed()){
+    os << var.value();
+  }else{
+    os << var.description();
+  }
+  return os;
+}

Linear_Programming/Variable.h

@@ -3,6 +3,7 @@
 
 #include <string>
 #include <iostream>
+#include <cassert>
 
 #include "../Common/integrality.h"
 
@@ -12,6 +13,8 @@ class Variable{
   integrality _integrality;
   std::pair<bool, double> _lower_bound;
   std::pair<bool, double> _upper_bound;
+
+  bool _fixed;
   double _value;
 public:
   Variable(const Variable& var);
@@ -23,9 +26,13 @@ public:
   integrality is_integral();
   std::pair<bool, double> lower_bound();
   std::pair<bool, double> upper_bound();
+
+  bool is_fixed();
   double& value();
 
   static bool VERBOSE_IDENT;
 };
 
+std::ostream& operator<<(std::ostream& os, Variable& var);
+
 #endif

Linear_Programming/lp_generator.h

@@ -23,7 +23,7 @@ using edge_constraint_data_generator = std::function<std::vector<constraint_data
 
 class lp_generator{
   template <typename N, typename E>
-  static std::map<Variable*, double> lhs_from_data(
+  static std::map<Variable*, Coefficient> lhs_from_data(
     size_t start_index, size_t length,
     lhs_constraint_data<N,E>& data,
     std::map<std::pair<Edge<N,E>*,E>, std::pair<Variable*, size_t>>& edge_variable_lookup,

Linear_Programming/lp_generator.ipp

 template <typename N, typename E>
-std::map<Variable*, double> lp_generator::lhs_from_data(size_t start_index, size_t length, lhs_constraint_data<N,E>& data, std::map<std::pair<Edge<N,E>*,E>, std::pair<Variable*, size_t>>& edge_variable_lookup, std::map<std::pair<Node<N,E>*,N>, std::pair<Variable*, size_t>>& node_variable_lookup){
-  std::map<Variable*, double> lhs;
+std::map<Variable*, Coefficient> lp_generator::lhs_from_data(size_t start_index, size_t length, lhs_constraint_data<N,E>& data, std::map<std::pair<Edge<N,E>*,E>, std::pair<Variable*, size_t>>& edge_variable_lookup, std::map<std::pair<Node<N,E>*,N>, std::pair<Variable*, size_t>>& node_variable_lookup){
+  std::map<Variable*, Coefficient> lhs;
 
   for(std::pair<std::pair<Node<N,E>*,N>, double> w : data.first){
     auto search_result = node_variable_lookup.find(w.first);
     if(search_result == node_variable_lookup.end()) throw std::range_error("lp_generator::lhs_from_data use of undeclared variable in constraint generation! (node)");
-    lhs.insert({search_result->second.first, w.second});
+    lhs.insert({search_result->second.first, {w.second}});
   }
   for(std::pair<std::pair<Edge<N,E>*,E>, double> w : data.second){
     auto search_result = edge_variable_lookup.find(w.first);
     if(search_result == edge_variable_lookup.end()) throw std::range_error("lp_generator::lhs_from_data use of undeclared variable in constraint generation! (edge)");
-    lhs.insert({search_result->second.first, w.second});
+    lhs.insert({search_result->second.first, {w.second}});
   }
 
   return lhs;

Makefile

@@ -54,7 +54,9 @@ LP := ./Linear_Programming/
 Common := ./Common_Types/
 Linear_Program_dep := Polyeder.o
 Polyeder_dep := Constraint.o
-Constraint_dep := Variable.o
+Constraint_dep := Variable.o Coefficient.o
+Coefficient_dep := Variable.o Monomial.o
+Monomial_dep := Variable.o
 maintenance_problem_generator_dep := Maintenance_Problem.o