diff --git a/examples/Gravity_test.cpp b/examples/Gravity_test.cpp index 4cd09779..fa1b30d3 100755 --- a/examples/Gravity_test.cpp +++ b/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 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"){ diff --git a/examples/Optimization/MINLP/LABS/LABS.cpp b/examples/Optimization/MINLP/LABS/LABS.cpp index dc502865..f3ffcf04 100644 --- a/examples/Optimization/MINLP/LABS/LABS.cpp +++ b/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> models; + for (int i = 0; i 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 s("s", -1, 1); var<> z("z", -1, 1); var 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 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 y_on_off("y_on_off"); // y_on_off = s - 2*y + 1; diff --git a/include/gravity/constant.h b/include/gravity/constant.h index a1c6a96e..d46d595f 100755 --- a/include/gravity/constant.h +++ b/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; } diff --git a/include/gravity/func.h b/include/gravity/func.h index 310f7a54..99b2d3bb 100755 --- a/include/gravity/func.h +++ b/include/gravity/func.h @@ -7945,6 +7945,13 @@ namespace gravity { shared_ptr p_new1; shared_ptr 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*)(&p1); + if(v->_range->first==-1 && v->_range->second==1){ + this->add_cst(constant(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 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>(p_it.first); + if(v->_range->first==-1 && v->_range->second==1){ + auto simplified = make_shared, 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; diff --git a/include/gravity/param.h b/include/gravity/param.h index f436a228..cb469531 100755 --- a/include/gravity/param.h +++ b/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; diff --git a/include/gravity/var.h b/include/gravity/var.h index 2a14156c..12847587 100644 --- a/include/gravity/var.h +++ b/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 res(*this); res._indices->_type = excl_; diff --git a/src/GurobiProgram.cpp b/src/GurobiProgram.cpp index dd95db93..8defdd33 100755 --- a/src/GurobiProgram.cpp +++ b/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);