summaryrefslogtreecommitdiff
path: root/omega/omega_lib/src/hull.cc
diff options
context:
space:
mode:
Diffstat (limited to 'omega/omega_lib/src/hull.cc')
-rw-r--r--omega/omega_lib/src/hull.cc1489
1 files changed, 1489 insertions, 0 deletions
diff --git a/omega/omega_lib/src/hull.cc b/omega/omega_lib/src/hull.cc
new file mode 100644
index 0000000..f1b0601
--- /dev/null
+++ b/omega/omega_lib/src/hull.cc
@@ -0,0 +1,1489 @@
+/*****************************************************************************
+ Copyright (C) 1994-2000 University of Maryland
+ Copyright (C) 2008 University of Southern California
+ Copyright (C) 2009-2010 University of Utah
+ All Rights Reserved.
+
+ Purpose:
+ Various hull calculations.
+
+ Notes:
+
+ History:
+ 06/15/09 ConvexRepresentation, Chun Chen
+ 11/25/09 RectHull, Chun Chen
+*****************************************************************************/
+
+#include <omega.h>
+#include <omega/farkas.h>
+#include <omega/hull.h>
+#include <basic/Bag.h>
+#include <basic/Map.h>
+#include <basic/omega_error.h>
+#include <list>
+#include <vector>
+#include <set>
+
+namespace omega {
+
+int hull_debug = 0;
+
+Relation ConvexHull(NOT_CONST Relation &R) {
+ Relation S = Approximate(consume_and_regurgitate(R));
+ if (!S.is_upper_bound_satisfiable())
+ return S;
+ if (S.has_single_conjunct())
+ return S;
+ return Farkas(Farkas(S,Basic_Farkas), Convex_Combination_Farkas);
+}
+
+Relation DecoupledConvexHull(NOT_CONST Relation &R) {
+ Relation S = Approximate(consume_and_regurgitate(R));
+ if (!S.is_upper_bound_satisfiable())
+ return S;
+ if (S.has_single_conjunct())
+ return S;
+ return Farkas(Farkas(S,Decoupled_Farkas), Convex_Combination_Farkas);
+}
+
+Relation AffineHull(NOT_CONST Relation &R) {
+ Relation S = Approximate(consume_and_regurgitate(R));
+ if (!S.is_upper_bound_satisfiable())
+ return S;
+ return Farkas(Farkas(S,Basic_Farkas), Affine_Combination_Farkas);
+}
+
+Relation LinearHull(NOT_CONST Relation &R) {
+ Relation S = Approximate(consume_and_regurgitate(R));
+ if (!S.is_upper_bound_satisfiable())
+ return S;
+ return Farkas(Farkas(S,Basic_Farkas), Linear_Combination_Farkas);
+}
+
+Relation ConicHull(NOT_CONST Relation &R) {
+ Relation S = Approximate(consume_and_regurgitate(R));
+ if (!S.is_upper_bound_satisfiable())
+ return S;
+ return Farkas(Farkas(S,Basic_Farkas), Positive_Combination_Farkas);
+}
+
+
+Relation FastTightHull(NOT_CONST Relation &input_R, NOT_CONST Relation &input_H) {
+ Relation R = Approximate(consume_and_regurgitate(input_R));
+ Relation H = Approximate(consume_and_regurgitate(input_H));
+
+ if (hull_debug) {
+ fprintf(DebugFile,"[ Computing FastTightHull of:\n");
+ R.prefix_print(DebugFile);
+ fprintf(DebugFile,"given known hull of:\n");
+ H.prefix_print(DebugFile);
+ }
+
+ if (!H.has_single_conjunct()) {
+ if (hull_debug)
+ fprintf(DebugFile, "] bailing out of FastTightHull, known hull not convex\n");
+ return H;
+ }
+
+ if (!H.is_obvious_tautology()) {
+ R = Gist(R,copy(H));
+ R.simplify(1,0);
+ }
+
+ if (R.has_single_conjunct()) {
+ R = Intersection(R,H);
+ if (hull_debug) {
+ fprintf(DebugFile, "] quick easy answer to FastTightHull\n");
+ R.prefix_print(DebugFile);
+ }
+ return R;
+ }
+ if (R.has_local(coefficient_of_constant_term)) {
+ if (hull_debug) {
+ fprintf(DebugFile, "] Can't handle recursive application of Farkas lemma\n");
+ }
+ return H;
+ }
+
+ if (hull_debug) {
+ fprintf(DebugFile,"Gist of R given H is:\n");
+ R.prefix_print(DebugFile);
+ }
+
+ if (1) {
+ Set<Variable_ID> vars;
+ int conjuncts = 0;
+ for (DNF_Iterator s(R.query_DNF()); s.live(); s.next()) {
+ conjuncts++;
+ for (Variable_ID_Iterator v(*((*s)->variables())); v.live(); v++) {
+ bool found = false;
+ for (EQ_Iterator eq = (*s)->EQs(); eq.live(); eq.next())
+ if ((*eq).get_coef(*v) != 0) {
+ if (!found) vars.insert(*v);
+ found = true;
+ break;
+ }
+ if (!found)
+ for (GEQ_Iterator geq = (*s)->GEQs(); geq.live(); geq.next())
+ if ((*geq).get_coef(*v) != 0) {
+ if (!found) vars.insert(*v);
+ found = true;
+ break;
+ }
+ }
+ }
+
+
+ // We now know which variables appear in R
+ if (hull_debug) {
+ fprintf(DebugFile,"Variables we need a better hull on are: ");
+ foreach(v,Variable_ID,vars,
+ fprintf(DebugFile," %s",v->char_name()));
+ fprintf(DebugFile,"\n");
+ }
+ Conjunct *c = H.single_conjunct();
+ int total=0;
+ int copied = 0;
+ for (EQ_Iterator eq = c->EQs(); eq.live(); eq.next()) {
+ total++;
+ foreach(v,Variable_ID,vars,
+ if ((*eq).get_coef(v) != 0) {
+ R.and_with_EQ(*eq);
+ copied++;
+ break; // out of variable loop
+ }
+ );
+ }
+ for (GEQ_Iterator geq = c->GEQs(); geq.live(); geq.next()) {
+ total++;
+ foreach(v,Variable_ID,vars,
+ if ((*geq).get_coef(v) != 0) {
+ R.and_with_GEQ(*geq);
+ copied++;
+ break; // out of variable loop
+ }
+ );
+ }
+ if (copied < total) {
+ R = Approximate(R);
+
+ if (hull_debug) {
+ fprintf(DebugFile,"Decomposed relation, copied only %d of %d constraints\n",copied,total);
+ fprintf(DebugFile,"Original R:\n");
+ R.prefix_print(DebugFile);
+ fprintf(DebugFile,"Known hull:\n");
+ H.prefix_print(DebugFile);
+ fprintf(DebugFile,"New R:\n");
+ R.prefix_print(DebugFile);
+ }
+ }
+ }
+
+ Relation F = Farkas(copy(R), Basic_Farkas, true);
+ if (hull_debug)
+ fprintf(DebugFile,"Farkas Difficulty = " coef_fmt "\n", farkasDifficulty);
+ if (farkasDifficulty > 260) {
+ if (hull_debug) {
+ fprintf(DebugFile, "] bailing out, farkas is way too complex\n");
+ fprintf(DebugFile,"Farkas:\n");
+ F.prefix_print(DebugFile);
+ }
+ return H;
+ }
+ else if (farkasDifficulty > 130) {
+ // Bail out
+ if (hull_debug) {
+ fprintf(DebugFile, coef_fmt " non-zeros in original farkas\n", farkasDifficulty);
+ }
+ Relation tmp = Farkas(R, Decoupled_Farkas, true);
+
+ if (hull_debug) {
+ fprintf(DebugFile, coef_fmt " non-zeros in decoupled farkas\n", farkasDifficulty);
+ }
+ if (farkasDifficulty > 260) {
+ if (hull_debug) {
+ fprintf(DebugFile, "] bailing out, farkas is way too complex\n");
+ fprintf(DebugFile,"Farkas:\n");
+ F.prefix_print(DebugFile);
+ }
+ return H;
+ }
+ else {
+ if (farkasDifficulty > 130)
+ R = Intersection(H, Farkas(tmp, Affine_Combination_Farkas, true));
+ else R = Intersection(H,
+ Intersection(Farkas(tmp, Convex_Combination_Farkas, true),
+ Farkas(F, Affine_Combination_Farkas, true)));
+ if (hull_debug) {
+ fprintf(DebugFile, "] bailing out, farkas is too complex, using affine hull\n");
+ fprintf(DebugFile,"Farkas:\n");
+ F.prefix_print(DebugFile);
+ fprintf(DebugFile,"Affine hull:\n");
+ R.prefix_print(DebugFile);
+ }
+ return R;
+ }
+ }
+
+ R = Intersection(H, Farkas(F, Convex_Combination_Farkas, true));
+ if (hull_debug) {
+ fprintf(DebugFile, "] Result of FastTightHull:\n");
+ R.prefix_print(DebugFile);
+ }
+ return R;
+}
+
+
+
+namespace {
+ bool parallel(const GEQ_Handle &g1, const GEQ_Handle &g2) {
+ for(Constr_Vars_Iter cvi(g1, false); cvi; cvi++) {
+ coef_t c1 = (*cvi).coef;
+ coef_t c2 = g2.get_coef((*cvi).var);
+ if (c1 != c2) return false;
+ }
+ {
+ for(Constr_Vars_Iter cvi(g2, false); cvi; cvi++) {
+ coef_t c1 = g1.get_coef((*cvi).var);
+ coef_t c2 = (*cvi).coef;
+ if (c1 != c2) return false;
+ }
+ }
+ return true;
+ }
+
+
+ bool hull(const EQ_Handle &e, const GEQ_Handle &g, coef_t &hull) {
+ int sign = 0;
+ for(Constr_Vars_Iter cvi(e, false); cvi; cvi++) {
+ coef_t c1 = (*cvi).coef;
+ coef_t c2 = g.get_coef((*cvi).var);
+ if (sign == 0) sign = (c1*c2>=0?1:-1);
+ if (sign*c1 != c2) return false;
+ }
+ assert(sign != 0);
+ {
+ for(Constr_Vars_Iter cvi(g, false); cvi; cvi++) {
+ coef_t c1 = e.get_coef((*cvi).var);
+ coef_t c2 = (*cvi).coef;
+ if (sign*c1 != c2) return false;
+ }
+ }
+ hull = max(sign * e.get_const(), g.get_const());
+ if (hull_debug) {
+ fprintf(DebugFile,"Hull of:\n %s\n", e.print_to_string().c_str());
+ fprintf(DebugFile," %s\n", g.print_to_string().c_str());
+ fprintf(DebugFile,"is " coef_fmt "\n\n",hull);
+ }
+ return true;
+ }
+
+ bool eq(const EQ_Handle &e1, const EQ_Handle &e2) {
+ int sign = 0;
+ for(Constr_Vars_Iter cvi(e1, false); cvi; cvi++) {
+ coef_t c1 = (*cvi).coef;
+ coef_t c2 = e2.get_coef((*cvi).var);
+ if (sign == 0) sign = (c1*c2>=0?1:-1);
+ if (sign*c1 != c2) return false;
+ }
+ assert(sign != 0);
+ {
+ for(Constr_Vars_Iter cvi(e2, false); cvi; cvi++) {
+ coef_t c1 = e1.get_coef((*cvi).var);
+ coef_t c2 = (*cvi).coef;
+ if (sign*c1 != c2) return false;
+ }
+ }
+ return sign * e1.get_const() == e2.get_const();
+ }
+}
+
+
+// This function is deprecated!!!
+Relation QuickHull(Relation &R) {
+ Tuple<Relation> Rs(1);
+ Rs[1] = R;
+ return QuickHull(Rs);
+}
+
+
+// This function is deprecated!!!
+Relation QuickHull(Tuple<Relation> &Rs) {
+ assert(!Rs.empty());
+
+ // if (Rs.size() == 1) return Rs[1];
+
+ Tuple<Relation> l_Rs;
+ for (int i = 1; i <= Rs.size(); i++)
+ for (DNF_Iterator c(Rs[i].query_DNF()); c; c++) {
+ Relation r = Relation(Rs[i], c.curr());
+ l_Rs.append(Approximate(r));
+ }
+
+ if (l_Rs.size() == 1)
+ return l_Rs[1];
+
+ Relation result = Relation::True(Rs[1]);
+ result.copy_names(Rs[1]);
+
+ use_ugly_names++;
+
+ if (hull_debug > 1)
+ for (int i = 1; i <= l_Rs.size(); i++) {
+ fprintf(DebugFile,"#%d \n",i);
+ l_Rs[i].prefix_print(DebugFile);
+ }
+
+
+// Relation R = copy(Rs[1]);
+// for (int i = 2; i <= Rs.size(); i++)
+// R = Union(R,copy(Rs[i]));
+
+// #if 0
+// if (!R.is_set()) {
+// if (R.n_inp() == R.n_out()) {
+// Relation AC = DeltasToRelation(Hull(Deltas(copy(R),
+// min(R.n_inp(),R.n_out()))),
+// R.n_inp(),R.n_out());
+// Relation dH = Hull(Domain(copy(R)),false);
+// Relation rH = Hull(Range(copy(R)),false);
+// result = Intersection(AC,Cross_Product(dH,rH));
+// }
+// else {
+// Relation dH = Hull(Domain(copy(R)),false);
+// Relation rH = Hull(Range(copy(R)),false);
+// result = Cross_Product(dH,rH);
+// assert(Must_Be_Subset(copy(R),copy(result)));
+// }
+// }
+
+// #endif
+
+ Conjunct *first;
+ l_Rs[1] = EQs_to_GEQs(l_Rs[1]);
+ first = l_Rs[1].single_conjunct();
+ for (GEQ_Iterator candidate(first->GEQs()); candidate.live(); candidate.next()) {
+ coef_t maxConstantTerm = (*candidate).get_const();
+ bool found = true;
+ if (hull_debug > 1) {
+ fprintf(DebugFile,"searching for bound on:\n %s\n", (*candidate).print_to_string().c_str());
+ }
+ for (int i = 2; i <= l_Rs.size(); i++) {
+ Conjunct *other = l_Rs[i].single_conjunct();
+ bool found_for_i = false;
+ for (GEQ_Iterator target(other->GEQs()); target.live(); target.next()) {
+ if (hull_debug > 2) {
+ fprintf(DebugFile,"candidate:\n %s\n", (*candidate).print_to_string().c_str());
+ fprintf(DebugFile,"target:\n %s\n", (*target).print_to_string().c_str());
+ }
+ if (parallel(*candidate,*target)) {
+ if (hull_debug > 1)
+ fprintf(DebugFile,"Found bound:\n %s\n", (*target).print_to_string().c_str());
+ maxConstantTerm = max(maxConstantTerm,(*target).get_const());
+ found_for_i = true;
+ break;
+ }
+ }
+ if (!found_for_i) {
+ for (EQ_Iterator target_e(other->EQs()); target_e.live(); target_e.next()) {
+ coef_t h;
+ if (hull(*target_e,*candidate,h)) {
+ if (hull_debug > 1)
+ fprintf(DebugFile,"Found bound of " coef_fmt ":\n %s\n", h, (*target_e).print_to_string().c_str());
+ maxConstantTerm = max(maxConstantTerm,h);
+ found_for_i = true;
+ break;
+ }
+ };
+ if (!found_for_i) {
+ if (hull_debug > 1) {
+ fprintf(DebugFile,"No bound found in:\n");
+ fprintf(DebugFile, "%s", l_Rs[i].print_with_subs_to_string().c_str());
+ }
+ //if nothing found
+ found = false;
+ break;
+ }
+ }
+ }
+
+ if (found) {
+ GEQ_Handle h = result.and_with_GEQ();
+ copy_constraint(h,*candidate);
+ if (hull_debug > 1)
+ fprintf(DebugFile,"Setting constant term to " coef_fmt " in\n %s\n", maxConstantTerm, h.print_to_string().c_str());
+ h.update_const(maxConstantTerm - (*candidate).get_const());
+ if (hull_debug > 1)
+ fprintf(DebugFile,"Updated constraint is\n %s\n", h.print_to_string().c_str());
+ }
+ }
+
+
+ for (EQ_Iterator candidate_eq(first->EQs()); candidate_eq.live(); candidate_eq.next()) {
+ bool found = true;
+ for (int i = 2; i <= l_Rs.size(); i++) {
+ Conjunct *C = l_Rs[i].single_conjunct();
+ bool found_for_i = false;
+
+ for (EQ_Iterator target(C->EQs()); target.live(); target.next()) {
+ if (eq(*candidate_eq,*target)) {
+ found_for_i = true;
+ break;
+ }
+ }
+ if (!found_for_i) {
+ //if nothing found
+ found = false;
+ break;
+ }
+ }
+
+ if (found) {
+ EQ_Handle h = result.and_with_EQ();
+ copy_constraint(h,*candidate_eq);
+ if (hull_debug > 1)
+ fprintf(DebugFile,"Adding eq constraint: %s\n", h.print_to_string().c_str());
+ }
+ }
+
+ use_ugly_names--;
+ if (hull_debug > 1) {
+ fprintf(DebugFile,"quick hull is of:");
+ result.print_with_subs(DebugFile);
+ }
+ return result;
+}
+
+
+// Relation Hull2(Tuple<Relation> &Rs, Tuple<int> &active) {
+// assert(Rs.size() == active.size() && Rs.size() > 0);
+
+// Tuple<Relation> l_Rs;
+// for (int i = 1; i <= Rs.size(); i++)
+// if (active[i])
+// l_Rs.append(copy(Rs[i]));
+
+// if (l_Rs.size() == 0)
+// return Relation::False(Rs[1]);
+
+// try {
+// Relation r = l_Rs[1];
+// for (int i = 2; i <= l_Rs.size(); i++) {
+// r = Union(r, copy(l_Rs[i]));
+// r.simplify();
+// }
+
+// // Relation F = Farkas(r, Basic_Farkas, true);
+// // if (farkasDifficulty >= 500)
+// // throw std::overflow_error("loop convex hull too complicated.");
+// // F = Farkas(F, Convex_Combination_Farkas, true);
+// return Farkas(Farkas(r, Basic_Farkas, true), Convex_Combination_Farkas, true);
+// }
+// catch (std::overflow_error) {
+// return QuickHull(l_Rs);
+// }
+// }
+
+
+namespace {
+ void printRs(Tuple<Relation> &Rs) {
+ fprintf(DebugFile,"Rs:\n");
+ for (int i = 1; i <= Rs.size(); i++)
+ fprintf(DebugFile,"#%d : %s\n",i,
+ Rs[i].print_with_subs_to_string().c_str());
+ }
+}
+
+Relation BetterHull(Tuple<Relation> &Rs, bool stridesAllowed, bool checkSubsets,
+ NOT_CONST Relation &input_knownHull = Relation::Null()) {
+ Relation knownHull = consume_and_regurgitate(input_knownHull);
+ static int OMEGA_WHINGE = -1;
+ if (OMEGA_WHINGE < 0) {
+ OMEGA_WHINGE = getenv("OMEGA_WHINGE") ? atoi(getenv("OMEGA_WHINGE")) : 0;
+ }
+ assert(!Rs.empty());
+ if (Rs.size() == 1) {
+ if (stridesAllowed) return Rs[1];
+ else return Approximate(Rs[1]);
+ }
+
+ if (checkSubsets) {
+ Tuple<bool> live(Rs.size());
+ if (hull_debug) {
+ fprintf(DebugFile,"Checking subsets in hull computation:\n");
+ printRs(Rs);
+ }
+ int i;
+ for(i=1;i <=Rs.size(); i++) live[i] = true;
+ for(i=1;i <=Rs.size(); i++)
+ for(int j=1;j <=Rs.size(); j++) if (i != j && live[j]) {
+ if (hull_debug) fprintf(DebugFile,"checking %d Is_Obvious_Subset %d\n",i,j);
+ if (Is_Obvious_Subset(copy(Rs[i]),copy(Rs[j]))) {
+ if (hull_debug) fprintf(DebugFile,"yes...\n");
+ live[i] = false;
+ break;
+ }
+ }
+ for(i=1;i <=Rs.size(); i++) if (!live[i]) {
+ if (i < Rs.size()) {
+ Rs[i] = Rs[Rs.size()];
+ live[i] = live[Rs.size()];
+ }
+ Rs[Rs.size()] = Relation();
+ Rs.delete_last();
+ i--;
+ }
+ }
+ Relation hull;
+ if (hull_debug) {
+ fprintf(DebugFile,"Better Hull:\n");
+ printRs(Rs);
+ fprintf(DebugFile,"known hull: %s\n", knownHull.print_with_subs_to_string().c_str());
+ }
+ if (knownHull.is_null()) hull = QuickHull(Rs);
+ else hull = Intersection(QuickHull(Rs),knownHull);
+ // for (int i = 1; i <= Rs.size(); i++)
+ // hull = RectHull(Union(hull, copy(Rs[i])));
+ // hull = Intersection(hull, knownHull);
+ hull.simplify();
+ if (hull_debug) {
+ fprintf(DebugFile,"quick hull: %s\n", hull.print_with_subs_to_string().c_str());
+ }
+
+ Relation orig = Relation::False(Rs[1]);
+ int i;
+ for (i = 1; i <= Rs.size(); i++)
+ orig = Union(orig,copy(Rs[i]));
+
+ orig.simplify();
+
+ for (i = 1; i <= Rs.size(); i++) {
+ if (!hull.is_obvious_tautology()) Rs[i] = Gist(Rs[i],copy(hull));
+ Rs[i].simplify();
+ if (Rs[i].is_obvious_tautology()) return hull;
+ if (Rs[i].has_single_conjunct()) {
+ Rs[i] = EQs_to_GEQs(Rs[i]);
+ if (hull_debug) {
+ fprintf(DebugFile,"Checking for hull constraints in:\n %s\n", Rs[i].print_with_subs_to_string().c_str());
+ }
+ Conjunct *c = Rs[i].single_conjunct();
+ for (GEQ_Iterator g(c->GEQs()); g.live(); g.next()) {
+ Relation tmp = Relation::True(Rs[i]);
+ tmp.and_with_GEQ(*g);
+ if (!Difference(copy(orig),tmp).is_upper_bound_satisfiable())
+ hull.and_with_GEQ(*g);
+ }
+ for (EQ_Iterator e(c->EQs()); e.live(); e.next()) {
+ Relation tmp = Relation::True(Rs[i]);
+ tmp.and_with_EQ(*e);
+ if (!Difference(copy(orig),tmp).is_upper_bound_satisfiable())
+ hull.and_with_EQ(*e);
+ }
+ }
+ }
+
+ hull = FastTightHull(orig,hull);
+ assert(hull.has_single_conjunct());
+
+ if (stridesAllowed) return hull;
+ else return Approximate(hull);
+
+}
+
+
+
+Relation Hull(NOT_CONST Relation &S,
+ bool stridesAllowed,
+ int effort,
+ NOT_CONST Relation &knownHull) {
+ Relation R = consume_and_regurgitate(S);
+ R.simplify(1,0);
+ if (!R.is_upper_bound_satisfiable()) return R;
+ Tuple<Relation> Rs;
+ for (DNF_Iterator c(R.query_DNF()); c.live(); ) {
+ Rs.append(Relation(R,c.curr()));
+ c.next();
+ }
+ if (effort == 1)
+ return BetterHull(Rs,stridesAllowed,false,knownHull);
+ else
+ return QuickHull(Rs);
+}
+
+
+
+Relation Hull(Tuple<Relation> &Rs,
+ Tuple<int> &validMask,
+ int effort,
+ bool stridesAllowed,
+ NOT_CONST Relation &knownHull) {
+ // Use relation of index i only when validMask[i] != 0
+ Tuple<Relation> Rs2;
+ for(int i = 1; i <= Rs.size(); i++) {
+ if (validMask[i]) {
+ Rs[i].simplify();
+ for (DNF_Iterator c(Rs[i].query_DNF()); c.live(); ) {
+ Rs2.append(Relation(Rs[i],c.curr()));
+ c.next();
+ }
+ }
+ }
+ assert(effort == 0 || effort == 1);
+ if (effort == 1)
+ return BetterHull(Rs2,stridesAllowed,true,knownHull);
+ else
+ return QuickHull(Rs2);
+}
+
+
+// This function is deprecated!!!
+Relation CheckForConvexRepresentation(NOT_CONST Relation &R_In) {
+ Relation R = consume_and_regurgitate(R_In);
+ Relation h = Hull(copy(R));
+ if (!Difference(copy(h),copy(R)).is_upper_bound_satisfiable())
+ return h;
+ else
+ return R;
+}
+
+// This function is deprecated!!!
+Relation CheckForConvexPairs(NOT_CONST Relation &S) {
+ Relation R = consume_and_regurgitate(S);
+ Relation hull = FastTightHull(copy(R),Relation::True(R));
+ R.simplify(1,0);
+ if (!R.is_upper_bound_satisfiable() || R.number_of_conjuncts() < 2) return R;
+ Tuple<Relation> Rs;
+ for (DNF_Iterator c(R.query_DNF()); c.live(); ) {
+ Rs.append(Relation(R,c.curr()));
+ c.next();
+ }
+
+ bool *dead = new bool[Rs.size()+1];
+ int i;
+ for(i = 1; i<=Rs.size();i++) dead[i] = false;
+
+ for(i = 1; i<=Rs.size();i++)
+ if (!dead[i])
+ for(int j = i+1; j<=Rs.size();j++) if (!dead[j]) {
+ if (hull_debug) {
+ fprintf(DebugFile,"Comparing #%d and %d\n",i,j);
+ }
+ Relation U = Union(copy(Rs[i]),copy(Rs[j]));
+ Relation H_ij = FastTightHull(copy(U),copy(hull));
+ if (!Difference(copy(H_ij),U).is_upper_bound_satisfiable()) {
+ Rs[i] = H_ij;
+ dead[j] = true;
+ if (hull_debug) {
+ fprintf(DebugFile,"Combined them\n");
+ }
+ }
+ }
+ i = 1;
+ while(i<=Rs.size() && dead[i]) i++;
+ assert(i<=Rs.size());
+ R = Rs[i];
+ i++;
+ for(; i<=Rs.size();i++)
+ if (!dead[i])
+ R = Union(R,Rs[i]);
+ delete []dead;
+ return R;
+}
+
+//
+// Supporting functions for ConvexRepresentation
+//
+namespace {
+struct Interval {
+ std::list<std::pair<Relation, Relation> >::iterator pos;
+ coef_t lb;
+ coef_t ub;
+ bool change;
+ coef_t modulo;
+ Interval(std::list<std::pair<Relation, Relation> >::iterator pos_, coef_t lb_, coef_t ub_):
+ pos(pos_), lb(lb_), ub(ub_) {}
+ friend bool operator<(const Interval &a, const Interval &b);
+};
+
+bool operator<(const Interval &a, const Interval &b) {
+ return a.lb < b.lb;
+}
+
+struct Modulo_Interval {
+ coef_t modulo;
+ coef_t start;
+ coef_t size;
+ Modulo_Interval(coef_t modulo_, coef_t start_, coef_t size_):
+ modulo(modulo_), start(start_), size(size_) {}
+ friend bool operator<(const Interval &a, const Interval &b);
+};
+
+bool operator<(const Modulo_Interval &a, const Modulo_Interval &b) {
+ if (a.modulo == b.modulo) {
+ if (a.start == b.start)
+ return a.size < b.size;
+ else
+ return a.start < b.start;
+ }
+ else
+ return a.modulo < b.modulo;
+}
+
+void merge_intervals(std::list<Interval> &intervals, coef_t modulo, std::list<std::pair<Relation, Relation> > &Rs, std::list<std::pair<Relation, Relation> >::iterator orig) {
+ // normalize intervals
+ for (std::list<Interval>::iterator i = intervals.begin(); i != intervals.end(); i++) {
+ (*i).modulo = modulo;
+ (*i).change = false;
+ if ((*i).ub - (*i).lb + 1>= modulo) {
+ (*i).lb = 0;
+ (*i).ub = modulo - 1;
+ }
+ else if ((*i).ub < 0 || (*i).lb >= modulo) {
+ coef_t range = (*i).ub - (*i).lb;
+ (*i).lb = int_mod((*i).lb, modulo);
+ (*i).ub = (*i).lb + range;
+ }
+ }
+
+ intervals.sort();
+
+ // merge neighboring intervals
+ std::list<Interval>::iterator p = intervals.begin();
+ while (p != intervals.end()) {
+ std::list<Interval>::iterator q = p;
+ q++;
+ while (q != intervals.end()) {
+ if ((*p).ub + 1 >= (*q).lb) {
+ Relation hull = ConvexHull(Union(copy((*(*p).pos).first), copy((*(*q).pos).first)));
+ Relation remainder = Difference(Difference(copy(hull), copy((*(*p).pos).first)), copy((*(*q).pos).first));
+ if (!remainder.is_upper_bound_satisfiable()) {
+ if ((*q).pos == orig)
+ std::swap((*p).pos, (*q).pos);
+ (*(*p).pos).first = hull;
+ (*p).ub = max((*p).ub, (*q).ub);
+ (*p).change = true;
+ Rs.erase((*q).pos);
+ q = intervals.erase(q);
+ }
+ else
+ break;
+ }
+ else
+ break;
+ }
+
+ bool p_moved = false;
+ q = p;
+ q++;
+ while (q != intervals.end()) {
+ if ((*q).ub >= modulo && int_mod((*q).ub, modulo) + 1 >= (*p).lb) {
+ Relation hull = ConvexHull(Union(copy((*(*p).pos).first), copy((*(*q).pos).first)));
+ Relation remainder = Difference(Difference(copy(hull), copy((*(*p).pos).first)), copy((*(*q).pos).first));
+ if (!remainder.is_upper_bound_satisfiable()) {
+ if ((*p).pos == orig)
+ std::swap((*p).pos, (*q).pos);
+ (*(*q).pos).first = hull;
+ coef_t t = (*p).ub - int_mod((*q).ub, modulo);
+ if (t > 0)
+ (*q).ub = (*q).ub + t;
+ (*q).change = true;
+ Rs.erase((*p).pos);
+ p = intervals.erase(p);
+ p_moved = true;
+ break;
+ }
+ else
+ q++;
+ }
+ else
+ q++;
+ }
+
+ if (!p_moved)
+ p++;
+ }
+
+ // merge by reducing the strengh of modulo
+ std::list<Modulo_Interval> modulo_intervals;
+ coef_t max_distance = modulo/2;
+ for (std::list<Interval>::iterator p = intervals.begin(); p != intervals.end(); p++) {
+ if ((*p).lb >= max_distance)
+ break;
+
+ coef_t size = (*p).ub - (*p).lb;
+
+ std::list<Interval>::iterator q = p;
+ q++;
+ while (q != intervals.end()) {
+ coef_t distance = (*q).lb - (*p).lb;
+ if (distance > max_distance)
+ break;
+
+ if ((*q).ub - (*q).lb != size || int_mod(modulo, distance) != 0) {
+ q++;
+ continue;
+ }
+
+ int num_reduced = 0;
+ coef_t looking_for = int_mod((*p).lb, distance);
+ for (std::list<Interval>::iterator k = intervals.begin(); k != intervals.end(); k++) {
+ if ((*k).lb == looking_for && (*k).ub - (*k).lb == size) {
+ num_reduced++;
+ looking_for += distance;
+ if (looking_for >= modulo)
+ break;
+ }
+ else if ((*k).lb <= looking_for && (*k).ub >= looking_for + size) {
+ looking_for += distance;
+ if (looking_for >= modulo)
+ break;
+ }
+ else if ((*k).lb > looking_for)
+ break;
+ }
+
+ if (looking_for >= modulo && num_reduced > 1)
+ modulo_intervals.push_back(Modulo_Interval(distance, int_mod((*p).lb, distance), size));
+
+ q++;
+ }
+ }
+
+ modulo_intervals.sort();
+
+ // remove redundant reduced-strength intervals
+ std::list<Modulo_Interval>::iterator p2 = modulo_intervals.begin();
+ while (p2 != modulo_intervals.end()) {
+ std::list<Modulo_Interval>::iterator q2 = p2;
+ q2++;
+ while (q2 != modulo_intervals.end()) {
+ if ((*p2).modulo == (*q2).modulo && (*p2).start == (*q2).start)
+ q2 = modulo_intervals.erase(q2);
+ else if (int_mod((*q2).modulo, (*p2).modulo) == 0 &&
+ (*p2).start == int_mod((*q2).start, (*p2).modulo) &&
+ (*p2).size >= (*q2).size)
+ q2 = modulo_intervals.erase(q2);
+ else
+ q2++;
+ }
+ p2++;
+ }
+
+ // replace original intervals with new reduced-strength ones
+ for (std::list<Modulo_Interval>::iterator i = modulo_intervals.begin(); i != modulo_intervals.end(); i++) {
+ std::vector<Relation *> candidates;
+ int num_replaced = 0;
+ for (std::list<Interval>::iterator j = intervals.begin(); j != intervals.end(); j++)
+ if (int_mod((*j).modulo, (*i).modulo) == 0 &&
+ (*j).ub - (*j).lb >= (*i).size &&
+ (int_mod((*j).lb, (*i).modulo) == (*i).start ||
+ int_mod((*j).ub, (*i).modulo) == (*i).start + (*i).size)) {
+ candidates.push_back(&((*(*j).pos).first));
+ if (int_mod((*j).lb, (*i).modulo) == (*i).start &&
+ (*j).ub - (*j).lb == (*i).size)
+ num_replaced++;
+ }
+ if (num_replaced <= 1)
+ continue;
+
+ Relation R = copy(*candidates[0]);
+ for (size_t k = 1; k < candidates.size(); k++)
+ R = Union(R, copy(*candidates[k]));
+ Relation hull = ConvexHull(copy(R));
+ Relation remainder = Difference(copy(hull), copy(R));
+ if (!remainder.is_upper_bound_satisfiable()) {
+ std::list<Interval>::iterator replaced_one = intervals.end();
+ for (std::list<Interval>::iterator j = intervals.begin(); j != intervals.end();)
+ if (int_mod((*j).modulo, (*i).modulo) == 0 &&
+ (*j).ub - (*j).lb >= (*i).size &&
+ (int_mod((*j).lb, (*i).modulo) == (*i).start ||
+ int_mod((*j).ub, (*i).modulo) == (*i).start + (*i).size)) {
+ if (int_mod((*j).lb, (*i).modulo) == (*i).start &&
+ (*j).ub - (*j).lb == (*i).size) {
+ if (replaced_one == intervals.end()) {
+ (*(*j).pos).first = hull;
+ (*j).lb = int_mod((*j).lb, (*i).modulo);
+ (*j).ub = int_mod((*j).ub, (*i).modulo);
+ (*j).modulo = (*i).modulo;
+ (*j).change = true;
+ replaced_one = j;
+ j++;
+ }
+ else {
+ if ((*j).pos == orig) {
+ std::swap((*replaced_one).pos, (*j).pos);
+ (*(*replaced_one).pos).first = (*(*j).pos).first;
+ }
+ Rs.erase((*j).pos);
+ j = intervals.erase(j);
+ }
+ }
+ else {
+ if (int_mod((*j).lb, (*i).modulo) == (*i).start)
+ (*j).lb = (*j).lb + (*i).size + 1;
+ else
+ (*j).ub = (*j).ub - (*i).size - 1;
+ (*j).change = true;
+ j++;
+ }
+ }
+ else
+ j++;
+ }
+ }
+}
+} // namespace
+
+
+//
+// Simplify a union of sets/relations to a minimal (may not be
+// optimal) number of convex regions. It intends to replace
+// CheckForConvexRepresentation and CheckForConvexPairs functions.
+//
+Relation ConvexRepresentation(NOT_CONST Relation &R) {
+ Relation l_R = copy(R);
+ if (!l_R.is_upper_bound_satisfiable() || l_R.number_of_conjuncts() < 2)
+ return R;
+
+ // separate each conjunct into smooth convex region and holes
+ std::list<std::pair<Relation, Relation> > Rs; // pair(smooth region, hole condition)
+ for (DNF_Iterator c(l_R.query_DNF()); c.live(); c++) {
+ Relation r1 = Relation(l_R, c.curr());
+ Relation r2 = Approximate(copy(r1));
+ r1 = Gist(r1, copy(r2));
+ Rs.push_back(std::make_pair(r2, r1));
+ }
+
+ try {
+ bool change = true;
+ while (change) {
+ change = false;
+
+ std::list<std::pair<Relation, Relation> >::iterator i = Rs.begin();
+ while (i != Rs.end()) {
+ // find regions with identical hole conditions to merge
+ {
+ std::list<std::pair<Relation, Relation> >::iterator j = i;
+ j++;
+ while (j != Rs.end()) {
+ if (!Difference(copy((*i).second), copy((*j).second)).is_upper_bound_satisfiable() &&
+ !Difference(copy((*j).second), copy((*i).second)).is_upper_bound_satisfiable()) {
+ if (Must_Be_Subset(copy((*j).first), copy((*i).first))) {
+ j = Rs.erase(j);
+ }
+ else if (Must_Be_Subset(copy((*i).first), copy((*j).first))) {
+ (*i).first = (*j).first;
+ j = Rs.erase(j);
+ change = true;
+ }
+ else {
+ Relation r;
+ bool already_use_recthull = false;
+ try {
+ // chun's debug
+ // throw std::runtime_error("dfdf");
+
+ r = ConvexHull(Union(copy((*i).first), copy((*j).first)));
+ }
+ catch (const std::overflow_error &e) {
+ r = RectHull(Union(copy((*i).first), copy((*j).first)));
+ already_use_recthull = true;
+ }
+ retry_recthull:
+ Relation r2 = Difference(Difference(copy(r), copy((*i).first)), copy((*j).first));
+ if (!r2.is_upper_bound_satisfiable()) { // convex hull is tight
+ (*i).first = r;
+ j = Rs.erase(j);
+ change = true;
+ }
+ else {
+ if (!already_use_recthull) {
+ r = RectHull(Union(copy((*i).first), copy((*j).first)));
+ already_use_recthull = true;
+ goto retry_recthull;
+ }
+ else
+ j++;
+ }
+ }
+ }
+ else
+ j++;
+ }
+ }
+
+ // find identical smooth regions as candidates for hole merge
+ std::list<std::list<std::pair<Relation, Relation> >::iterator> s;
+ for (std::list<std::pair<Relation, Relation> >::iterator j = Rs.begin(); j != Rs.end(); j++)
+ if (j != i) {
+ if (!Intersection(Difference(copy((*i).first), copy((*j).first)), copy((*j).second)).is_upper_bound_satisfiable() &&
+ !Intersection(Difference(copy((*j).first), copy((*i).first)), copy((*i).second)).is_upper_bound_satisfiable())
+ s.push_back(j);
+ }
+
+ if (s.size() != 0) {
+ // convert hole condition c1*x1+c2*x2+... = c*alpha+d to a pair of inequalities
+ (*i).second = EQs_to_GEQs((*i).second, false);
+
+ // find potential wildcards that can be used for hole conditions
+ std::set<Variable_ID> nonsingle_wild;
+ for (EQ_Iterator ei((*i).second.single_conjunct()); ei; ei++)
+ if ((*ei).has_wildcards())
+ for (Constr_Vars_Iter cvi(*ei, true); cvi; cvi++)
+ nonsingle_wild.insert(cvi.curr_var());
+ for (GEQ_Iterator gei((*i).second.single_conjunct()); gei; gei++)
+ if ((*gei).has_wildcards()) {
+ Constr_Vars_Iter cvi(*gei, true);
+ Constr_Vars_Iter cvi2 = cvi;
+ cvi2++;
+ if (cvi2) {
+ nonsingle_wild.insert(cvi.curr_var());
+ for (; cvi2; cvi2++)
+ nonsingle_wild.insert(cvi2.curr_var());
+ }
+ }
+
+ // find hole condition in c*alpha+d1<=c1*x1+c2*x2+...<=c*alpha+d2 format
+ for (GEQ_Iterator gei((*i).second.single_conjunct()); gei; gei++)
+ if ((*gei).has_wildcards()) {
+ coef_t c;
+ Variable_ID v;
+ {
+ Constr_Vars_Iter cvi(*gei, true);
+ v = cvi.curr_var();
+ c = cvi.curr_coef();
+ if (c < 0 || nonsingle_wild.find(v) != nonsingle_wild.end())
+ continue;
+ }
+
+ coef_t lb = posInfinity;
+ for (GEQ_Iterator gei2((*i).second.single_conjunct()); gei2; gei2++) {
+ if (!(*gei2 == *gei) && (*gei2).get_coef(v) != 0) {
+ if (lb != posInfinity) {
+ nonsingle_wild.insert(v);
+ break;
+ }
+
+ bool match = true;
+ for (Constr_Vars_Iter cvi2(*gei); cvi2; cvi2++)
+ if (cvi2.curr_coef() != -((*gei2).get_coef(cvi2.curr_var()))) {
+ match = false;
+ break;
+ }
+ if (match)
+ for (Constr_Vars_Iter cvi2(*gei2); cvi2; cvi2++)
+ if (cvi2.curr_coef() != -((*gei).get_coef(cvi2.curr_var()))) {
+ match = false;
+ break;
+ }
+ if (!match) {
+ nonsingle_wild.insert(v);
+ break;
+ }
+
+ lb = -(*gei2).get_const();
+ }
+ }
+
+ if (nonsingle_wild.find(v) != nonsingle_wild.end())
+ continue;
+
+ Relation stride_cond = Relation::True((*i).second);
+ F_Exists *f_exists = stride_cond.and_with_and()->add_exists();
+ Variable_ID e = f_exists->declare();
+ F_And *f_root = f_exists->add_and();
+ GEQ_Handle h1 = f_root->add_GEQ();
+ GEQ_Handle h2 = f_root->add_GEQ();
+ for (Constr_Vars_Iter cvi2(*gei); cvi2; cvi2++) {
+ Variable_ID v = cvi2.curr_var();
+ switch (v->kind()) {
+ case Wildcard_Var:
+ h1.update_coef(e, cvi2.curr_coef());
+ h2.update_coef(e, -cvi2.curr_coef());
+ break;
+ case Global_Var: {
+ Global_Var_ID g = v->get_global_var();
+ Variable_ID v2;
+ if (g->arity() == 0)
+ v2 = stride_cond.get_local(g);
+ else
+ v2 = stride_cond.get_local(g, v->function_of());
+ h1.update_coef(v2, cvi2.curr_coef());
+ h2.update_coef(v2, -cvi2.curr_coef());
+ break;
+ }
+ default:
+ h1.update_coef(v, cvi2.curr_coef());
+ h2.update_coef(v, -cvi2.curr_coef());
+ }
+ }
+ h1.update_const((*gei).get_const());
+ h2.update_const(-lb);
+
+ stride_cond.simplify();
+ Relation other_cond = Gist(copy((*i).second), copy(stride_cond));
+
+ // find regions with potential mergeable stride condition with this one
+ std::list<Interval> intervals;
+ intervals.push_back(Interval(i, lb, (*gei).get_const()));
+
+ for (std::list<std::list<std::pair<Relation, Relation> >::iterator>::iterator j = s.begin(); j != s.end(); j++)
+ if (Must_Be_Subset(copy((**j).second), copy(other_cond))) {
+ Relation stride_cond2 = Gist(copy((**j).second), copy(other_cond));
+
+ // interval can be removed
+ if (stride_cond2.is_obvious_tautology()) {
+ intervals.push_back(Interval(*j, 0, c-1));
+ continue;
+ }
+
+ stride_cond2 = EQs_to_GEQs(stride_cond2, false);
+ coef_t lb, ub;
+ GEQ_Iterator gei2(stride_cond2.single_conjunct());
+ coef_t sign = 0;
+ for (Constr_Vars_Iter cvi(*gei2, true); cvi; cvi++)
+ if (sign != 0) {
+ sign = 0;
+ break;
+ }
+ else if (cvi.curr_coef() == c)
+ sign = 1;
+ else if (cvi.curr_coef() == -c)
+ sign = -1;
+ else {
+ sign = 0;
+ break;
+ }
+ if (sign == 0)
+ continue;
+
+ bool match = true;
+ for (Constr_Vars_Iter cvi(*gei2); cvi; cvi++) {
+ Variable_ID v = cvi.curr_var();
+ if (v->kind() == Wildcard_Var)
+ continue;
+ else if (v->kind() == Global_Var) {
+ Global_Var_ID g = v->get_global_var();
+ if (g->arity() == 0)
+ v = (*i).second.get_local(g);
+ else
+ v = (*i).second.get_local(g, v->function_of());
+ }
+
+ if (cvi.curr_coef() != sign * (*gei).get_coef(v)) {
+ match = false;
+ break;
+ }
+ }
+ if (!match)
+ continue;
+
+ for (Constr_Vars_Iter cvi(*gei); cvi; cvi++) {
+ Variable_ID v = cvi.curr_var();
+ if (v->kind() == Wildcard_Var)
+ continue;
+ else if (v->kind() == Global_Var) {
+ Global_Var_ID g = v->get_global_var();
+ if (g->arity() == 0)
+ v = stride_cond2.get_local(g);
+ else
+ v = stride_cond2.get_local(g, v->function_of());
+ }
+
+ if (cvi.curr_coef() != sign * (*gei2).get_coef(v)) {
+ match = false;
+ break;
+ }
+ }
+ if (!match)
+ continue;
+ if (sign > 0)
+ ub = (*gei2).get_const();
+ else
+ lb = -(*gei2).get_const();
+
+ gei2++;
+ if (!gei2)
+ continue;
+
+ coef_t sign2 = 0;
+ for (Constr_Vars_Iter cvi(*gei2, true); cvi; cvi++)
+ if (sign2 != 0) {
+ sign2 = 0;
+ break;
+ }
+ else if (cvi.curr_coef() == c)
+ sign2 = 1;
+ else if (cvi.curr_coef() == -c)
+ sign2 = -1;
+ else {
+ sign2 = 0;
+ break;
+ }
+ if (sign2 != -sign)
+ continue;
+
+ for (Constr_Vars_Iter cvi(*gei2); cvi; cvi++) {
+ Variable_ID v = cvi.curr_var();
+ if (v->kind() == Wildcard_Var)
+ continue;
+ else if (v->kind() == Global_Var) {
+ Global_Var_ID g = v->get_global_var();
+ if (g->arity() == 0)
+ v = (*i).second.get_local(g);
+ else
+ v = (*i).second.get_local(g, v->function_of());
+ }
+
+ if (cvi.curr_coef() != sign2 * (*gei).get_coef(v)) {
+ match = false;
+ break;
+ }
+ }
+ if (!match)
+ continue;
+
+ for (Constr_Vars_Iter cvi(*gei); cvi; cvi++) {
+ Variable_ID v = cvi.curr_var();
+ if (v->kind() == Wildcard_Var)
+ continue;
+ else if (v->kind() == Global_Var) {
+ Global_Var_ID g = v->get_global_var();
+ if (g->arity() == 0)
+ v = stride_cond2.get_local(g);
+ else
+ v = stride_cond2.get_local(g, v->function_of());
+ }
+
+ if (cvi.curr_coef() != sign2 * (*gei2).get_coef(v)) {
+ match = false;
+ break;
+ }
+ }
+ if (!match)
+ continue;
+ if (sign2 > 0)
+ ub = (*gei2).get_const();
+ else
+ lb = -(*gei2).get_const();
+
+ gei2++;
+ if (gei2)
+ continue;
+
+ intervals.push_back(Interval(*j, lb, ub));
+ }
+
+ merge_intervals(intervals, c, Rs, i);
+
+ // make current region the last one being updated
+ bool invalid = false;
+ for (std::list<Interval>::iterator ii = intervals.begin(); ii != intervals.end(); ii++)
+ if ((*ii).change && (*ii).pos == i) {
+ invalid = true;
+ intervals.push_back(*ii);
+ intervals.erase(ii);
+ break;
+ }
+
+ // update hole condition for each region
+ for (std::list<Interval>::iterator ii = intervals.begin(); ii != intervals.end(); ii++)
+ if ((*ii).change) {
+ change = true;
+
+ if ((*ii).ub - (*ii).lb + 1 >= (*ii).modulo)
+ (*(*ii).pos).second = copy(other_cond);
+ else {
+ Relation stride_cond = Relation::True((*i).second);
+ F_Exists *f_exists = stride_cond.and_with_and()->add_exists();
+ Variable_ID e = f_exists->declare();
+ F_And *f_root = f_exists->add_and();
+ GEQ_Handle h1 = f_root->add_GEQ();
+ GEQ_Handle h2 = f_root->add_GEQ();
+ for (Constr_Vars_Iter cvi2(*gei); cvi2; cvi2++) {
+ Variable_ID v = cvi2.curr_var();
+ switch (v->kind()) {
+ case Wildcard_Var:
+ h1.update_coef(e, (*ii).modulo);
+ h2.update_coef(e, -(*ii).modulo);
+ break;
+ case Global_Var: {
+ Global_Var_ID g = v->get_global_var();
+ Variable_ID v2;
+ if (g->arity() == 0)
+ v2 = stride_cond.get_local(g);
+ else
+ v2 = stride_cond.get_local(g, v->function_of());
+ h1.update_coef(v2, cvi2.curr_coef());
+ h2.update_coef(v2, -cvi2.curr_coef());
+ break;
+ }
+ default:
+ h1.update_coef(v, cvi2.curr_coef());
+ h2.update_coef(v, -cvi2.curr_coef());
+ }
+ }
+ h1.update_const((*ii).ub);
+ h2.update_const(-(*ii).lb);
+
+ (*(*ii).pos).second = Intersection(copy(other_cond), stride_cond);
+ (*(*ii).pos).second.simplify();
+ }
+ }
+
+ if (invalid)
+ break;
+ }
+ }
+ i++;
+ }
+ }
+ }
+ catch (const presburger_error &e) {
+ throw e;
+ }
+
+ Relation R2 = Relation::False(l_R);
+ for (std::list<std::pair<Relation, Relation> >::iterator i = Rs.begin(); i != Rs.end(); i++)
+ R2 = Union(R2, Intersection((*i).first, (*i).second));
+ R2.simplify(0, 1);
+
+ return R2;
+}
+
+
+//
+// Use gist and value range to calculate a quick rectangular hull. It
+// intends to replace all hull calculations (QuickHull, BetterHull,
+// FastTightHull) beyond the method of ConvexHull (dual
+// representations). In the future, it will support max(...)-like
+// upper bound. So RectHull complements ConvexHull in two ways: first
+// for relations that ConvexHull gets too complicated, second for
+// relations where different conjuncts have different symbolic upper
+// bounds.
+//
+Relation RectHull(NOT_CONST Relation &Rel) {
+ Relation R = Approximate(consume_and_regurgitate(Rel));
+ if (!R.is_upper_bound_satisfiable())
+ return R;
+ if (R.has_single_conjunct())
+ return R;
+
+ std::vector<std::string> input_names(R.n_inp());
+ for (int i = 1; i <= R.n_inp(); i++)
+ input_names[i-1] = R.input_var(i)->name();
+ std::vector<std::string> output_names(R.n_out());
+ for (int i = 1; i <= R.n_out(); i++)
+ output_names[i-1] = R.output_var(i)->name();
+
+ DNF_Iterator c(R.query_DNF());
+ Relation r = Relation(R, c.curr());
+ c++;
+ std::vector<std::pair<coef_t, coef_t> > bounds1(R.n_inp());
+ std::vector<std::pair<coef_t, coef_t> > bounds2(R.n_out());
+ {
+ Relation t = Project_Sym(copy(r));
+ t.simplify();
+ for (int i = 1; i <= R.n_inp(); i++) {
+ Tuple<Variable_ID> v;
+ for (int j = 1; j <= R.n_inp(); j++)
+ if (j != i)
+ v.append(r.input_var(j));
+ for (int j = 1; j <= R.n_out(); j++)
+ v.append(r.output_var(j));
+ Relation t2 = Project(copy(t), v);
+ t2.query_variable_bounds(t2.input_var(i), bounds1[i-1].first, bounds1[i-1].second);
+ }
+ for (int i = 1; i <= R.n_out(); i++) {
+ Tuple<Variable_ID> v;
+ for (int j = 1; j <= R.n_out(); j++)
+ if (j != i)
+ v.append(r.output_var(j));
+ for (int j = 1; j <= R.n_inp(); j++)
+ v.append(r.input_var(j));
+ Relation t2 = Project(copy(t), v);
+ t2.query_variable_bounds(t2.output_var(i), bounds2[i-1].first, bounds2[i-1].second);
+ }
+ }
+
+ while (c.live()) {
+ Relation r2 = Relation(R, c.curr());
+ c++;
+ Relation x = Gist(copy(r), Gist(copy(r), copy(r2), 1), 1);
+ if (Difference(copy(r2), copy(x)).is_upper_bound_satisfiable())
+ x = Relation::True(R);
+ Relation y = Gist(copy(r2), Gist(copy(r2), copy(r), 1), 1);
+ if (Difference(copy(r), copy(y)).is_upper_bound_satisfiable())
+ y = Relation::True(R);
+ r = Intersection(x, y);
+
+ {
+ Relation t = Project_Sym(copy(r2));
+ t.simplify();
+ for (int i = 1; i <= R.n_inp(); i++) {
+ Tuple<Variable_ID> v;
+ for (int j = 1; j <= R.n_inp(); j++)
+ if (j != i)
+ v.append(r2.input_var(j));
+ for (int j = 1; j <= R.n_out(); j++)
+ v.append(r2.output_var(j));
+ Relation t2 = Project(copy(t), v);
+ coef_t lbound, ubound;
+ t2.query_variable_bounds(t2.input_var(i), lbound, ubound);
+ bounds1[i-1].first = min(bounds1[i-1].first, lbound);
+ bounds1[i-1].second = max(bounds1[i-1].second, ubound);
+ }
+ for (int i = 1; i <= R.n_out(); i++) {
+ Tuple<Variable_ID> v;
+ for (int j = 1; j <= R.n_out(); j++)
+ if (j != i)
+ v.append(r2.output_var(j));
+ for (int j = 1; j <= R.n_inp(); j++)
+ v.append(r2.input_var(j));
+ Relation t2 = Project(copy(t), v);
+ coef_t lbound, ubound;
+ t2.query_variable_bounds(t2.output_var(i), lbound, ubound);
+ bounds2[i-1].first = min(bounds2[i-1].first, lbound);
+ bounds2[i-1].second = max(bounds2[i-1].second, ubound);
+ }
+ }
+
+ Relation r3(R.n_inp(), R.n_out());
+ F_And *f_root = r3.add_and();
+ for (int i = 1; i <= R.n_inp(); i++) {
+ if (bounds1[i-1].first != -posInfinity) {
+ GEQ_Handle h = f_root->add_GEQ();
+ h.update_coef(r3.input_var(i), 1);
+ h.update_const(-bounds1[i-1].first);
+ }
+ if (bounds1[i-1].second != posInfinity) {
+ GEQ_Handle h = f_root->add_GEQ();
+ h.update_coef(r3.input_var(i), -1);
+ h.update_const(bounds1[i-1].second);
+ }
+ }
+ for (int i = 1; i <= R.n_out(); i++) {
+ if (bounds2[i-1].first != -posInfinity) {
+ GEQ_Handle h = f_root->add_GEQ();
+ h.update_coef(r3.output_var(i), 1);
+ h.update_const(-bounds2[i-1].first);
+ }
+ if (bounds2[i-1].second != posInfinity) {
+ GEQ_Handle h = f_root->add_GEQ();
+ h.update_coef(r3.output_var(i), -1);
+ h.update_const(bounds2[i-1].second);
+ }
+ }
+ r = Intersection(r, r3);
+ r.simplify();
+ }
+
+ for (int i = 1; i <= r.n_inp(); i++)
+ r.name_input_var(i, input_names[i-1]);
+ for (int i = 1; i <= r.n_out(); i++)
+ r.name_output_var(i, output_names[i-1]);
+ r.setup_names();
+ return r;
+}
+
+} // namespace
generated by cgit v1.2.3-70-g09d2 (git 2.48.1) at 2025-02-02 23:19:56 +0000