/*****************************************************************************
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