working in simplification in polynomes

examples/Gravity_test.cpp

@@ -22,9 +22,24 @@
 using namespace std;
 using namespace gravity;
 
-TEST_CASE("testing MISDP solvers"){
-    var<> x1("x1", -10, 150);
-
+TEST_CASE("testing simplification of polynomial terms with -1,1 integers squared"){
+    var<int> s("s", -1, 1);
+    s.in(range(0,3));
+    s.exclude_zero();
+
+    auto f1 = s[0]*s[0];
+    f1.print();
+    CHECK(f1.is_constant());
+    CHECK(f1==1);
+    auto f2 = s[1]*s[2]*s[1];
+    f2.print();
+    CHECK(f2.is_linear());
+    auto f3 = s[0]*s[1]*s[2]*s[1];
+    f3.print();
+    CHECK(f3.is_quadratic());
+    auto f4 = s[0]*s[0]*s[2]*s[1];
+    f4.print();
+    CHECK(f4.is_quadratic());
 }
 
 TEST_CASE("testing range update"){

examples/Optimization/MINLP/LABS/LABS.cpp

@@ -21,28 +21,43 @@ int main(int argc, char * argv[]){
     if(argc>=2){
         n=stoi(argv[1]);
     }
-    DebugOn("Optimizing with N = " << n << endl);
+    DebugOn("Optimizing for N = " << n << endl);
+    bool new_algorithm = false;
+    if(new_algorithm){
+        vector<Model<>> models;
+        for (int i = 0; i<n-1; i++) {
+            Model<> M("LABS_"+to_string(n)+"_"+to_string(i));
+
+            models.push_back(M);
+        }
+        return 0;
+    }
+
+
+
     Model<> M_obj("M_obj_LABS_"+to_string(n));
     Model<> M("LABS_"+to_string(n));
-    var<> s("s", -1, 1);
+    var<int> s("s", -1, 1);
     var<> z("z", -1, 1);
     var<int> y("y", 0, 1);
     var<> cs("cs", pos_);
     var<> c("c");
     indices s_ids = range(0,n-1);
-
+    indices c_ids = range(1,n-1);
     bool unconstrained = false;
     if(unconstrained){
-        M_obj.add(y.in(s_ids));
-        func<> obj;        
+        M_obj.add(s.in(s_ids), y.in(s_ids));
+        M.add(c.in(c_ids));
         for (int k = 1; k<=n-1; k++) {
             func<> cterm;
             for (int i = 0; i<n-k; i++) {
-                cterm += (2*y(to_string(i))-1)*(2*y(to_string(i+k))-1);
+                cterm += (s(to_string(i)))*(s(to_string(i+k)));
             }
-            obj += pow(cterm,2);
+            Constraint<> C_sq_k("C_sq_k");
+            C_sq_k = c(k) - pow(cterm,2);
+            M_obj.add(C_sq_k.in(c_ids)>=0);
         }
-        M_obj.min(obj);
+        M_obj.min(sum(c));
         M_obj.print();
         solver<> g_sol(M_obj,ipopt);
         g_sol.run();
@@ -51,7 +66,6 @@ int main(int argc, char * argv[]){
         M_obj.print_solution();
     }
     else{
-        indices c_ids = range(1,n-1);
         indices z_ids("z_ids");
         string idx;
         for (int i = 0; i<n-1; i++) {
@@ -62,7 +76,8 @@ int main(int argc, char * argv[]){
         }
         M.add(s.in(s_ids), c.in(c_ids), z.in(z_ids));
         M.add(y.in(s_ids));
-
+//        s.set_lb("0",1);
+//        y.set_lb("0",1);
 
     //    Constraint<> y_on_off("y_on_off");
     //    y_on_off = s - 2*y + 1;

include/gravity/constant.h

@@ -71,6 +71,7 @@ namespace gravity {
 
         bool                            _is_transposed = false; /**< True if the constant is transposed */
         bool                            _is_vector = false; /**< True if the constant is a vector or matrix */
+        bool                            _exclude_zero = false; /**< True if the constant cannot assume a zero value */
         size_t                          _dim[2] = {1,1}; /*< dimension of current object */
 
         bool                            _polar = false; /**< True in case this is a complex number with a polar representation, rectangular representation if false */
@@ -371,7 +372,8 @@ namespace gravity {
         constant& operator=(const constant& c) {
             _type = c._type;
             _is_transposed = c._is_transposed;
-            _is_vector = c._is_vector;                        
+            _is_vector = c._is_vector;
+            _exclude_zero = c._exclude_zero;
             _val = c._val;
             return *this;
         }
@@ -381,6 +383,7 @@ namespace gravity {
             update_type();
             _is_transposed = c._is_transposed;
             _is_vector = c._is_vector;
+            _exclude_zero = c._exclude_zero;
             _val = c._val;
             return *this;
         }

include/gravity/func.h

@@ -7945,6 +7945,13 @@ namespace gravity {
             shared_ptr<param_> p_new1;
             shared_ptr<param_> p_new2;
             _evaluated=false;
+            if(p1.is_var() && p2.is_var() && p1.is_integer() &&  p2.is_integer() && p1._exclude_zero && ps1==ps2){/* If this is p1^2 and p1 is a -1,1 binary, simplify expression */
+                auto v = (var<int>*)(&p1);
+                if(v->_range->first==-1 && v->_range->second==1){
+                    this->add_cst(constant<type>(1));
+                    return false;
+                }
+            }
             if (_ftype <= lin_ && p1.is_var()) {
                 _ftype = quad_;
             }
@@ -8127,6 +8134,7 @@ namespace gravity {
             auto pair_it = _pterms->find(name);
             auto p = l.begin()->first;
             shared_ptr<param_> pnew;
+            auto initial_ftype =_ftype;
             if (_ftype <= quad_ && p->is_var()) {
                 _ftype = pol_;
             }
@@ -8159,6 +8167,32 @@ namespace gravity {
                     newv = true;
                     for (auto& p_it:*newl) {
                         if (p_it.first->get_name(false,false)==s) {
+                            if(p_it.first->is_var() && p_it.first->is_integer() &&  p_it.first->_exclude_zero && p_it.second%2==1){/* If this is p^2 and p is a -1,1 binary, simplify expression */
+                                auto v = static_pointer_cast<var<int>>(p_it.first);
+                                if(v->_range->first==-1 && v->_range->second==1){
+                                    auto simplified = make_shared<list<std::pair<shared_ptr<param_>, int>>>();
+                                    for (auto& p_it2:*newl) {
+                                        if(p_it2!=p_it)
+                                            simplified->push_back(p_it2);
+                                    }
+                                    if(simplified->size()==1){
+                                        auto pt = simplified->front();
+                                        if(pt.second==1){
+                                            _ftype = initial_ftype;
+                                            this->insert(sign, coef, *pt.first);
+                                            return false;
+                                        }
+                                        if(pt.second==2){
+                                            _ftype = initial_ftype;
+                                            this->insert(sign, coef, *pt.first, *pt.first);
+                                            return false;
+                                        }
+                                    }
+                                    newl = simplified;
+                                    newv = false;
+                                    break;
+                                }
+                            }
                             p_it.second++;
                             newv = false;
                             break;

include/gravity/param.h

@@ -84,6 +84,7 @@ namespace gravity {
             _name = p._name;
             _is_transposed = p._is_transposed;
             _is_vector = p._is_vector;
+            _exclude_zero = p._exclude_zero;
             _is_angle = p._is_angle;
             _is_sqrmag = p._is_sqrmag;
             _is_conjugate = p._is_conjugate;
@@ -732,6 +733,7 @@ namespace gravity {
             res._name = _name;
             res._is_transposed = _is_transposed;
             res._is_vector = _is_vector;
+            res._exclude_zero = _exclude_zero;
             res._new = _new;
             res._is_relaxed = _is_relaxed;
             res._is_angle = _is_angle;
@@ -763,6 +765,7 @@ namespace gravity {
             _name = p._name;
             _is_transposed = p._is_transposed;
             _is_vector = p._is_vector;
+            _exclude_zero = p._exclude_zero;
             _new = p._new;
             _is_relaxed = p._is_relaxed;
             _is_angle = p._is_angle;
@@ -806,6 +809,7 @@ namespace gravity {
             _name = p._name;
             _is_transposed = p._is_transposed;
             _is_vector = p._is_vector;
+            _exclude_zero = p._exclude_zero;
             _new = p._new;
             _is_relaxed = p._is_relaxed;
             _is_angle = p._is_angle;
@@ -842,6 +846,7 @@ namespace gravity {
             _name = p._name;
             _is_transposed = p._is_transposed;
             _is_vector = p._is_vector;
+            _exclude_zero = p._exclude_zero;
             _new = p._new;
             _is_relaxed = p._is_relaxed;
             _is_angle = p._is_angle;

include/gravity/var.h

@@ -579,7 +579,9 @@ namespace gravity {
             return res;
         }
 
-
+        void exclude_zero() {
+            this->_exclude_zero = true;
+        }
         var excl(size_t index) {
             var<type> res(*this);
             res._indices->_type = excl_;

src/GurobiProgram.cpp

@@ -518,7 +518,7 @@ bool GurobiProgram::solve(bool relax, double mipgap, double time_limit){
     grb_mod->set(GRB_IntParam_OutputFlag,1);
 //    grb_mod->set(GRB_DoubleParam_Cutoff,59);
     //  grb_mod->set(GRB_IntParam_MinRelNodes,0);
-    //    grb_mod->set(GRB_DoubleParam_Heuristics, 1);
+//        grb_mod->set(GRB_DoubleParam_Heuristics,  0);
     //    grb_mod->set(GRB_DoubleParam_NoRelHeurTime, 5);
     //    grb_mod->set(GRB_DoubleParam_NoRelHeurWork, 5);
     //    grb_mod->set(GRB_IntParam_CutPasses,10000);