/*****************************************************************************
Copyright (C) 1994-2000 the Omega Project Team
Copyright (C) 2005-2011 Chun Chen
Copyright (C) 2009-2011 West Pomeranian University of Technology, Szczecin
All Rights Reserved.
Purpose:
All calculations of closure are now here.
Notes:
Related paper:
- "Transitive closure of infinite graphs and its applications",
Wayne Kelly, William Pugh, Evan Rosser and Tatiana Shpeisman, IJPP 1996.
- "Computing the Transitive Closure of a Union of Affine Integer Tuple
Relations", Anna Beletska, Denis Barthou, Wlodzimierz Bielecki and
Albert Cohen, COCOA 2009.
- "An Iterative Algorithm of Computing the Transitive Closure of a Union
of Parameterized Affine Integer Tuple Relations", Bielecki Wlodzimierz,
Klimek Tomasz, Palkowski Marek and Anna Beletska, COCOA 2010.
History:
12/27/09 move ConicClosure here, Chun Chen
01/19/11 new closure algorithms, Klimek Tomzsz
02/02/11 move VennDiagramFrom here, Chun Chen
*****************************************************************************/
#include <typeinfo>
#include <assert.h>
#include <omega.h>
#include <omega/hull.h>
#include <basic/Iterator.h>
#include <basic/List.h>
#include <basic/SimpleList.h>
namespace omega {
void InvestigateClosure(Relation r, Relation r_closure, Relation bounds);
void print_given_bounds(const Relation & R1, NOT_CONST Relation& input_Bounds);
#define printConjunctClosure (closure_presburger_debug & 0x1)
#define detailedClosureDebug (closure_presburger_debug & 0x2)
#ifdef TC_STATS
extern int clock_diff();
extern void start_clock();
FILE *statsfile;
int singles, totals=0;
#endif
int closure_presburger_debug = 0;
Relation VennDiagramForm(NOT_CONST Relation &Context_In,
Tuple<Relation> &Rs,
int next,
bool anyPositives,
int weight) {
Relation Context = consume_and_regurgitate(Context_In);
if (hull_debug) {
fprintf(DebugFile,"[VennDiagramForm, next = %d, anyPositives = %d, weight = %d \n", next,anyPositives,weight);
fprintf(DebugFile,"context:\n");
Context.prefix_print(DebugFile);
}
if (anyPositives && weight > 3) {
Context.simplify();
if (!Context.is_upper_bound_satisfiable()) {
if (hull_debug)
fprintf(DebugFile,"] not satisfiable\n");
return Context;
}
weight = 0;
}
if (next > Rs.size()) {
if (!anyPositives) {
if (hull_debug)
fprintf(DebugFile,"] no positives\n");
return Relation::False(Context);
}
Context.simplify();
if (hull_debug) {
fprintf(DebugFile,"] answer is:\n");
Context.prefix_print(DebugFile);
}
return Context;
}
Relation Pos = VennDiagramForm(Intersection(copy(Context),copy(Rs[next])),
Rs,
next+1,
true,
weight+2);
Relation Neg = VennDiagramForm(Difference(Context,copy(Rs[next])),
Rs,
next+1,
anyPositives,
weight+1);
if (hull_debug) {
fprintf(DebugFile,"] VennDiagramForm\n");
fprintf(DebugFile,"pos part:\n");
Pos.prefix_print(DebugFile);
fprintf(DebugFile,"neg part:\n");
Neg.prefix_print(DebugFile);
}
return Union(Pos,Neg);
}
Relation VennDiagramForm(Tuple<Relation> &Rs, NOT_CONST Relation &Context_In) {
Relation Context = consume_and_regurgitate(Context_In);
if (Context.is_null()) Context = Relation::True(Rs[1]);
if (hull_debug) {
fprintf(DebugFile,"Starting computation of VennDiagramForm\n");
fprintf(DebugFile,"Context:\n");
Context.prefix_print(DebugFile);
for(int i = 1; i <= Rs.size(); i++) {
fprintf(DebugFile,"#%d:\n",i);
Rs[i].prefix_print(DebugFile);
}
}
return VennDiagramForm(Context,Rs,1,false,0);
}
Relation VennDiagramForm(NOT_CONST Relation &R_In, NOT_CONST Relation &Context_In) {
Relation R = consume_and_regurgitate(R_In);
Relation Context = consume_and_regurgitate(Context_In);
Tuple<Relation> Rs;
for (DNF_Iterator c(R.query_DNF()); c.live(); ) {
Rs.append(Relation(R,c.curr()));
c.next();
}
return VennDiagramForm(Rs,Context);
}
Relation ConicClosure (NOT_CONST Relation &R) {
int n = R.n_inp();
if (n != R.n_out())
throw std::invalid_argument("conic closure must have the same input arity and output arity");
return DeltasToRelation(ConicHull(Deltas(R)), n, n);
}
bool is_lex_forward(Relation R) {
if(R.n_inp() != R.n_out()) {
fprintf(stderr, "relation has wrong inputs/outpts\n");
exit(1);
}
Relation forw(R.n_inp(), R.n_out());
F_Or * o = forw.add_or();
for(int a = 1; a <= forw.n_inp(); a++) {
F_And * andd = o->add_and();
GEQ_Handle g = andd->add_GEQ();
g.update_coef(input_var(a), -1);
g.update_coef(output_var(a), 1);
g.update_const(1);
for(int b = 1; b < a; b++) {
EQ_Handle e = andd->add_EQ();
e.update_coef(input_var(a),1);
e.update_coef(output_var(a),-1);
}
}
Relation test = Difference(R, forw);
return !test.is_upper_bound_satisfiable();
}
static Relation compose_n(NOT_CONST Relation &input_r, int n) {
Relation r = consume_and_regurgitate(input_r);
if (n == 1)
return r;
else
return Composition(r, compose_n(copy(r), n-1));
} /* compose_n */
Relation approx_closure(NOT_CONST Relation &input_r, int n) {
Relation r = consume_and_regurgitate(input_r);
Relation r_closure;
r_closure=r;
int i;
for(i=2; i<=n; i++)
r_closure=Union(r_closure,compose_n(copy(r), n));
r_closure = Union(r_closure, Relation::Unknown(r_closure));
return r_closure;
} /* approx_closure */
static bool is_closure_itself(NOT_CONST Relation &r) {
return Must_Be_Subset(Composition(copy(r),copy(r)),copy(r));
}
/*****
* get a D form of the Relation (single conjunct).
* D = {[ i_1,i_2,...,i_m] -> [j_1, j_2, ..., j_m ] :
* (forall p, 1<= p <= m) L_p <= j_p - i_p <= U_p &&
* j_p - i_p == M_p alpha_p};
* Right now only wildcards that are in stride constraints are treated.
*****/
Relation get_D_form (Relation & R) {
Relation D(R.n_inp(), R.n_out());
R.make_level_carried_to(R.n_inp());
assert(R.has_single_conjunct());
int n_zero=0;
for (DNF_Iterator d(R.query_DNF()); d.live(); d.next())
n_zero=d.curr()->query_guaranteed_leading_0s();
Relation Diff=Deltas(copy(R));
if (detailedClosureDebug) {
fprintf(DebugFile, "The relation projected onto differencies is:\n");
Diff.print_with_subs(DebugFile);
}
/* now form D */
int i;
coef_t l,u;
F_And * N = D.add_and();
GEQ_Handle g;
for (i=1; i<=Diff.n_set(); i++) {
Diff.query_variable_bounds(Diff.set_var(i), l,u);
/* if (i== n_zero+1 && l==negInfinity)
l=1; */
if (l!=negInfinity) {
g=N->add_GEQ();
g.update_coef(D.input_var(i),-1);
g.update_coef(D.output_var(i),1);
g.update_const(-l);
g.finalize();
}
if (u!=posInfinity) {
g=N->add_GEQ();
g.update_coef(D.input_var(i),1);
g.update_coef(D.output_var(i),-1);
g.update_const(u);
g.finalize();
}
}
/* add all stride constrains if they do exist */
Conjunct *c = Diff.single_conjunct();
if (c->locals().size()>0) {// there are local variables
// now go through all the equalities
coef_t coef=0;
int pos=0;
for (EQ_Iterator eq = c->EQs(); eq.live(); eq.next()) {
// constraint is in stride form if it has 2 vars,
// one of which is wildcard. Count number if vars and wildcard vars
int nwild=0,nvar=0;
for (Constr_Vars_Iter cvi(*eq, false); cvi; cvi++) {
if ((*cvi).var->kind() == Global_Var)
continue;
else if ((*cvi).var->kind() == Wildcard_Var) {
coef=(*cvi).coef;
nwild++;
}
else
pos=(*cvi).var->get_position();
nvar++;
}
if (nvar==2 && nwild==1) { //stride constraint
EQ_Handle e=N->add_stride(coef);
e.update_coef(D.input_var(pos),-1);
e.update_coef(D.output_var(pos),1);
e.finalize();
}
}
} // end search of stride constrains
D.finalize();
D.simplify();
return D;
} /* end get_D_form */
/****
* get relation A x A describing a region of domain and range:
* A=Hull(Domain(R), Range(R)) intersection IterationSpace
* returns cross product A x A
***/
Relation form_region(const Relation &R, const Relation& IterationSpace) {
Relation H=Union(Domain(copy(R)), Range(copy(R)));
H.simplify(1,1);
H = EQs_to_GEQs(H);
H=Hull(H);
Relation A=Intersection(H, copy(IterationSpace));
Relation A1=A;
return Cross_Product(A,A1);
}
Relation form_region1(const Relation &R, const Relation& IterationSpace) {
Relation Dom=Intersection(Domain(copy(R)), copy(IterationSpace));
Relation Ran=Intersection(Range(copy(R)), copy(IterationSpace));
return Cross_Product(Dom,Ran);
}
/****
* Check if we can use D instead of R
* i.e. D intersection (A cross A) is Must_Be_Subset of R
***/
bool isD_OK(Relation &R, Relation &D, Relation &AxA) {
Relation B=Intersection(copy(D), copy(AxA));
B.simplify();
if (detailedClosureDebug) {
fprintf(DebugFile, "Intersection of D and AxA is:\n");
B.print_with_subs(DebugFile);
}
assert (Must_Be_Subset(copy(R),copy(B)));
return Must_Be_Subset(B, copy(R));
}
/****
* check if the constraint is a stride one. Here we say that an equality
* constraint is a stride constraint if it has exatly one wildcard.
* The function returns number of the wildcards in the constraint.
* So if we know that constraint is from the relation in D form, then
* it cannot have more than 1 wildcard variables, and the result of
* this functions can be treated as bool.
***/
static int is_stride(const EQ_Handle &eq) {
int n=0;
for (Constr_Vars_Iter cvi(eq,true); cvi; cvi++)
n++;
return n;
}
/*****
* check if the constraint is in the form i_k' - i_k comp_op c
* return v - the number of the var and the type of the comp_op:
* 1 - >, -1 - <, 0 - not in the right form
* if this is equality constraint in the right form any 1 or -1 can be
* returned
******/
static coef_t is_constraint_in_D_form(Relation &r, const Constraint_Handle &h, int &v) {
v=-1;
coef_t c_out = 0;
for (int i = 1; i <= r.n_inp(); i++) {
coef_t c_in = h.get_coef(r.input_var(i));
if (c_in) {
if (v!=-1)
return 0;
v=i;
c_out = h.get_coef(r.output_var(i));
// special case for modular constraint -- by chun 04/02/2009
if (h.has_wildcards() && typeid(h) == typeid(EQ_Handle)) {
coef_t g = 0;
for (Constr_Vars_Iter cvi(h, true); cvi; cvi++)
g = gcd(g, abs(cvi.curr_coef()));
c_in = int_mod_hat(c_in, g);
c_out = int_mod_hat(c_out, g);
if (g == 2) {
if (c_in * c_out == 1) {
c_out = -1;
}
else
return 0;
}
else if (c_in * c_out != -1)
return 0;
}
// other cases
else if (c_in * c_out != -1)
return 0;
}
}
return c_out;
}
/***
* Check if relation is in the D form
* D = {[ i_1,i_2,...,i_m] -> [j_1, j_2, ..., j_m ] :
* (forall p, 1<= p <= m) L_p <= j_p - i_p <= U_p &&
* j_p - i_p == M_p alpha_p};
* Right now we do not check for multiple stride constraints for one var.
* Probably they cannot exist in simplified conjunct
* This function will be used in assertions
*****/
bool is_in_D_form(Relation & D) {
/* check that D has one conjunct */
if (! D.has_single_conjunct())
return false;
Conjunct * c=D.single_conjunct();
if (D.global_decls()->size() != 0) // there are symbolic vars
return false;
if (D.n_inp() != D.n_out())
return false;
int n=D.n_inp();
Tuple<int> bl(n), bu(n);
for (int i=1; i<= n; i++)
bl[i]=bu[i]=0;
int v;
coef_t res;
for (EQ_Iterator eq = c->EQs(); eq.live(); eq.next()) {
if ((res=is_constraint_in_D_form(D,*eq,v))==0)
return false;
int n_wild=is_stride(*eq);
if (n_wild>=2)
return false;
if (n_wild==0) { // not stride constraint
if (bl[v] || bu[v])
return false;
bl[v]=bu[v]=1;
}
}
for (GEQ_Iterator geq = c->GEQs(); geq.live(); geq.next()) {
if ((res=is_constraint_in_D_form(D,*geq,v))==0)
return false;
if ((res>0 && bl[v]) || (res<0 && bu[v]))
return false;
if (res>0)
bl[v]=1;
else
bu[v]=1;
}
return true;
}
#define get_D_plus_form(R) (get_D_closure(R,1))
#define get_D_star_form(R) (get_D_closure(R,0))
/****
* Get D+ or D* from the relation that is in D form
* To get D+ calculate:
* D+= {[i1, i2 .. i_m] -> {j1, j2, ..., j_m]:
* exists s s.t. s>=1 and
* (forall p, 1<= p <= m) L_p * s<= j_p - i_p <= U_p*s &&
* j_p - i_p == M_p alpha_p};
* To get D* calculate almost the same relation but s>=0.
* Parameter n is 1 for getting D+ and 0 for D*
****/
Relation get_D_closure(Relation & D, int n) {
assert (is_in_D_form(D));
assert(n==0 || n==1);
Conjunct *c=D.single_conjunct();
Relation R(D.n_inp(), D.n_out());
F_Exists * ex = R.add_exists();
Variable_ID s = ex->declare("s");
F_And * N = ex->add_and();
/* add s>=1 or s>=0 */
GEQ_Handle geq= N->add_GEQ();
geq.update_coef(s,1);
geq.update_const(-n);
geq.finalize();
/* copy and modify all the EQs */
for (EQ_Iterator j= c->EQs(); j.live(); j.next()) {
EQ_Handle eq=N->add_EQ();
copy_constraint(eq, *j);
// if it's stride constraint do not change it
if (!is_stride(*j)) {
/* eq is j_k -i_k = c, replace c buy s*c */
eq.update_coef(s, (*j).get_const());
eq.update_const(-(*j).get_const());
}
eq.finalize();
}
/* copy and modify all the GEQs */
for (GEQ_Iterator gi= c->GEQs(); gi.live(); gi.next()) {
geq=N->add_GEQ();
copy_constraint(geq, *gi);
/* geq is j_k -i_k >=c or i_k-j_k >=c, replace c buy s*c */
geq.update_coef(s,(*gi).get_const());
geq.update_const(-(*gi).get_const());
geq.finalize();
}
R.finalize();
if (detailedClosureDebug) {
fprintf(DebugFile, "Simplified D%c is:\n", n==1?'+':'*');
R.print_with_subs(DebugFile);
}
return R;
}
/***
* Check if we can easily calculate the D* (D* will be convex).
* We can calculate D* if all differences have both lower and upper
* bounds to be non -/+ infinity
***/
bool can_get_D_star_form(Relation &D) {
assert(is_in_D_form(D));
Conjunct *c=D.single_conjunct();
int n=D.n_inp();
Tuple<int> bl(n), bu(n);
int i;
for (i=1; i<=n; i++)
bl[i]=bu[i]=0;
for (EQ_Iterator eq = c->EQs(); eq.live(); eq.next()) {
// do not check stride constraints
if (!is_stride(*eq)) {
for (i=1; i<=n; i++) {
if ((*eq).get_coef(D.input_var(i)) !=0 )
bl[i]=bu[i]=1;
}
}
}
for (GEQ_Iterator geq = c->GEQs(); geq.live(); geq.next()) {
for (i=1; i<=n; i++) {
coef_t k;
if ((k=(*geq).get_coef(D.input_var(i))) != 0) {
if (k>0)
bu[i]=1;
else
bl[i]=1;
}
}
}
for (i=1; i<=n; i++)
if (!bl[i] || !bu[i])
return false;
return true;
}
/*****
* Check whether the relation intersect with identity or not
****/
bool does_intersect_with_identity(Relation &R) {
assert (R.n_inp() == R.n_out());
Relation I=Identity(R.n_inp());
Relation C=Intersection(I, copy(R));
return C.is_upper_bound_satisfiable();
}
bool does_include_identity(Relation &R) {
Relation I=Identity(R.n_inp());
return Must_Be_Subset(I, copy(R));
}
/*****
* Bill's closure: check if it is possible to calculate transitive closure
* of the relation using the Bill's algorithm.
* Return the transitive closure relation if it is possible and null relation
* otherwise
****/
bool Bill_closure(Relation &R, Relation& IterationSpace, Relation & R_plus, Relation & R_star) {
#ifdef TC_STATS
fprintf(statsfile,"start bill closure\n");
#endif
if (does_include_identity(R))
return false;
if (detailedClosureDebug) {
fprintf(DebugFile, "\nApplying Bill's method to calculate transitive closure\n");
}
// get D and AxA
Relation D=get_D_form(R);
if (detailedClosureDebug) {
fprintf(DebugFile,"\n D form for the relation:\n");
D.print_with_subs(DebugFile);
}
Relation AxA=form_region1(R, IterationSpace);
if (detailedClosureDebug) {
fprintf(DebugFile, "\n AxA for the relation:\n");
AxA.print_with_subs(DebugFile);
}
// compute R_+
R_plus=Intersection(get_D_plus_form(D), copy(AxA));
if (detailedClosureDebug) {
fprintf(DebugFile, "\nR_+= D+ intersection AxA is:\n");
R_plus.print_with_subs(DebugFile);
}
// compute R_*
R_star=Intersection(get_D_star_form(D), form_region(R,IterationSpace));
if (detailedClosureDebug) {
fprintf(DebugFile, "\nR_*= D* intersection AxA is:\n");
R_star.print_with_subs(DebugFile);
}
/* Check that R_+ is acyclic.
Given the way we constructed R_+, R_+=(R_+)+.
As a result it's enough to verify that R_+ intersection I = 0,
to prove that R_+ is acyclic.
*/
if (does_intersect_with_identity(R_plus)) {
if (detailedClosureDebug) {
fprintf(DebugFile,"R_+ is not acyclic.\n");
}
return false;
}
//Check R_+ - R is Must_Be_Subset of R o R_+
if (!Must_Be_Subset(Difference(copy(R_plus), copy(R)), Composition(copy(R), copy(R_plus)))) {
#if defined(TC_STATS)
fprintf(statsfile, "R_+ -R is not a Must_Be_Subset of R o R_+\n");
fprintf(statsfile, "Bill Method is not applicable\n");
#endif
return false;
}
if (detailedClosureDebug) {
fprintf(DebugFile, "R_+ -R is a Must_Be_Subset of R o R_+ - good\n");
}
// if we are here than all tests worked, and R_+ is transitive closure
// of R.
#if defined(TC_STATS)
fprintf(statsfile,"\nAll three tests succeeded -- exact closure found\n");
fprintf(statsfile, "Transitive closure is R_+\n");
#endif
// assert(isD_OK(R,D,AxA));
return true;
}
/**********************************************************************
* print the relation given the bounds on the iteration space
* If the bounds are unknown (Bounds is Null), then just print relation
* itself
****/
void print_given_bounds( const Relation& R1, NOT_CONST Relation& input_Bounds) {
Relation & Bounds = (Relation &)input_Bounds;
Relation r;
if (Bounds.is_null())
r=R1;
else
r = Gist(copy(R1),copy(Bounds),1);
r.print_with_subs(DebugFile);
}
/**********************************************************************
* Investigate closure:
* checks if the copmuted approximation on the Transitive closure
* is upper and lower bound. If it's both - it's exact.
* This function doesn't return any value. It's just prints a lot
* of debug output
* INPUT:
* r - relation
* r_closure - approximation on r+.
* F - iteration space
**********************************************************************/
void InvestigateClosure(Relation r, Relation r_closure, Relation F) {
Relation r3;
bool LB_res, UB_res;
if (!F.is_null())
F=Cross_Product(copy(F),copy(F));
fprintf(DebugFile, "\n\n--->investigating the closure of the relation:\n");
print_given_bounds(r,F);
fprintf(DebugFile, "\nComputed closure is:\n");
print_given_bounds(r_closure,F);
r3=Composition(copy(r),copy(r_closure));
r3.simplify(1,1);
r3=Union(r3,Composition(copy(r_closure),copy(r)));
r3.simplify(1,1);
r3=Union(r3,copy(r));
r3.simplify(1,1);
Relation remainder = Difference(copy(r3),copy(r_closure));
if (!F.is_null()) {
r3=Gist(r3,F,1);
}
r3.simplify(1,1);
if (!F.is_null()) {
r_closure=Gist(r_closure,F,1);
}
r_closure.simplify(1,1);
LB_res= Must_Be_Subset(copy(r_closure),copy(r3));
UB_res=Must_Be_Subset(copy(r3),copy(r_closure));
fprintf(DebugFile,"\nThe results of checking closure (gist) are:\n");
fprintf(DebugFile,"LB - %s, UB - %s\n", LB_res?"YES":"NO", UB_res?"YES":"NO");
if (!UB_res) {
remainder.simplify(2,2);
fprintf(DebugFile,"Dependences not included include:\n");
print_given_bounds(remainder,F);
}
}
/****
* Transitive closure of the relation containing single conjunct
****/
bool ConjunctTransitiveClosure (NOT_CONST Relation & input_R, Relation & IterationSpace, Relation & R_plus, Relation & R_star) {
Relation R = consume_and_regurgitate(input_R);
assert(R.has_single_conjunct());
if (printConjunctClosure) {
fprintf(DebugFile,"\nTaking closure of the single conjunct: [\n");
R.print_with_subs(DebugFile);
}
#ifdef TC_STATS
fprintf(statsfile,"start conjuncttransitiveclosure\n");
singles++;
#endif
if (is_closure_itself(copy(R))) {
#ifdef TC_STATS
fprintf(statsfile, "Relation is closure itself\n");
#endif
int ndim_all, ndim_domain;
R.dimensions(ndim_all,ndim_domain);
if (ndim_all == ndim_domain +1) {
Relation ispace = Cross_Product(Domain(copy(R)),Range(copy(R)));
Relation R_zero = Intersection(copy(ispace),Identity(R.n_inp()));
R_star = Hull(Union(copy(R),R_zero),true,1,ispace);
R_plus=R;
if (printConjunctClosure) {
fprintf(DebugFile, "\n] For this relation R+=R\n");
fprintf(DebugFile,"R*:\n");
R_star.print_with_subs(DebugFile);
}
return true;
}
else {
R_star=R;
R_plus=R;
if (printConjunctClosure) {
fprintf(DebugFile, "\n] For this relation R+=R, not appropriate for R*\n");
}
return false;
}
}
else {
bool done=false;
if (!IterationSpace.is_null()) {
// Bill's closure requires the information about Iteration Space.
// So if IterationSpace is NULL, i.e. unknown( e.g. when calling from parser,
// we do not do Bill's closure
done = Bill_closure(R, IterationSpace, R_plus, R_star);
#ifdef TC_STATS
fprintf(statsfile,"Bill closure is %sapplicable\n",done?"":"not ");
#endif
if (printConjunctClosure) {
if (!done)
fprintf(DebugFile, "Bill's closure is not applicable\n");
else {
fprintf(DebugFile, "Bill's closure is applicable\n");
fprintf (DebugFile, " For R:\n");
R.print_with_subs(DebugFile);
fprintf(DebugFile, "R+ is:\n");
R_plus.print_with_subs(DebugFile);
fprintf(DebugFile, "R* is:\n");
R_star.print_with_subs(DebugFile);
fprintf(DebugFile, "\n");
InvestigateClosure(R, R_plus, IterationSpace);
}
}
}
if (done) {
if (printConjunctClosure) {
fprintf(DebugFile, "]\n");
}
return true;
}
else {
// do and check approximate closure (several compositions)
R_plus = approx_closure(copy(R), 2);
#ifdef TC_STATS
fprintf(statsfile,"Approximating closure with 2 compositions\n");
#endif
if (printConjunctClosure) {
fprintf(DebugFile, "Doing approximate closure\n");
InvestigateClosure(R, R_plus, IterationSpace);
}
} //end else (!done after Bill Closure or Iteration space is NULL)
if (printConjunctClosure) {
fprintf(DebugFile, "]\n");
}
}
return false;
}
/*********************************************************************
* try to get conjunct transitive closure.
* it we can get it easy get it, return true.
* if not - return false
********************************************************************/
bool TryConjunctTransitiveClosure (NOT_CONST Relation & input_R, Relation & IterationSpace, Relation & R_plus) {
Relation R = consume_and_regurgitate(input_R);
assert(R.has_single_conjunct());
#ifdef TC_STATS
fprintf(statsfile,"start tryconjuncttransitiveclosure\n");
singles++;
#endif
if (printConjunctClosure) {
fprintf(DebugFile,"\nTrying to take closure of the single conjunct: [\n");
R.print_with_subs(DebugFile);
}
if (is_closure_itself(copy(R))) {
#ifdef TC_STATS
fprintf(statsfile, "Relation is closure itself, leave alone (try)\n");
#endif
if (printConjunctClosure)
fprintf(DebugFile, "\n ]The relation is closure itself. Leave it alone\n");
return false;
}
else {
bool done;
assert(!IterationSpace.is_null());
Relation R_star;
done = Bill_closure(R, IterationSpace, R_plus, R_star);
#ifdef TC_STATS
fprintf(statsfile, "Bill closure is %sapplicable (try)\n", done?"":"NOT ");
#endif
if (printConjunctClosure) {
if (!done)
fprintf(DebugFile, "]Bill's closure is not applicable\n");
else {
fprintf(DebugFile, "]Bill's closure is applicable\n");
fprintf (DebugFile, " For R:\n");
R.print_with_subs(DebugFile);
fprintf(DebugFile, "R+ is:\n");
R_plus.print_with_subs(DebugFile);
fprintf(DebugFile, "R* is:\n");
R_star.print_with_subs(DebugFile);
fprintf(DebugFile, "\n");
InvestigateClosure(R, R_plus, IterationSpace);
}
}
return done;
}
//return false;
}
bool Equal (const Relation & r1, const Relation & r2) {
bool res=Must_Be_Subset (copy(r1), copy(r2));
if (!res)
return false;
return Must_Be_Subset (copy(r2),copy(r1));
}
void appendClausesToList(Simple_List<Relation> &L, Relation &R) {
R.make_level_carried_to(R.n_inp());
R.simplify(2,2);
for(int depth = R.n_inp(); depth >= -1; depth--)
for (DNF_Iterator d(R.query_DNF()); d.live(); d.next())
if (d.curr()->query_guaranteed_leading_0s() == depth) {
L.append(Relation(R, d.curr()));
}
}
void printRelationList(Simple_List<Relation> &L) {
for (Simple_List_Iterator<Relation> li(L); li.live(); li.next()) {
li.curr().print_with_subs(DebugFile);
}
}
/****
* Transitive closure of the relation containing multiple conjuncts
* New (Bill's) version
***/
Relation TransitiveClosure0(NOT_CONST Relation &input_r, int maxExpansion, NOT_CONST Relation & input_IterationSpace) {
Relation r = consume_and_regurgitate(input_r);
Relation IterationSpace = consume_and_regurgitate(input_IterationSpace);
if (closure_presburger_debug)
fprintf(DebugFile, "\n\n[Transitive closure\n\n");
Relation result;
#ifdef TC_STATS
#define TC_RUNS 1
int in_conj = copy(r).query_DNF()->length();
totals++;
fprintf(statsfile,"%d closure run\n", totals);
if(is_in_D_form(copy(r)))
fprintf(statsfile, "Relation initially in D form\n");
else
fprintf(statsfile, "Relation initially NOT in D form\n");
if(is_lex_forward(copy(r)))
fprintf(statsfile, "Relation is initially lex forw\n");
else
fprintf(statsfile, "Relation is NOT initially lex forw\n");
start_clock();
for(int tc_loop = 1; tc_loop <= TC_RUNS; tc_loop++) {
singles = 0;
#endif
assert(!r.is_null());
assert(r.n_inp() == r.n_out());
if (r.max_ufs_arity() > 0) {
assert(r.max_ufs_arity() == 0 && "Can't take transitive closure with UFS yet.");
fprintf(stderr, "Can't take transitive closure with UFS yet.");
exit(1);
}
r.simplify(2,2);
if (!r.is_upper_bound_satisfiable()) {
#ifdef TC_STATS
int totalTime = clock_diff();
fprintf(statsfile, "Relation is unsatisfiable\n");
fprintf(statsfile, "input conj: %d output conj: %d #singe conj closures: %d time: %d\n",
in_conj, copy(result).query_DNF()->length(),
singles,
totalTime/TC_RUNS);
#endif
if (closure_presburger_debug)
fprintf(DebugFile, "]TC : relation is false\n");
return r;
}
IterationSpace = Hull(Union(Domain(copy(r)),Range(copy(r))), true, 1, IterationSpace);
if (detailedClosureDebug) {
fprintf(DebugFile, "r is:\n");
r.print_with_subs(DebugFile);
fprintf(DebugFile, "IS is:\n");
IterationSpace.print_with_subs(DebugFile);
}
Relation dom = Domain(copy(r));
dom.simplify(2,1);
Relation rng = Range(copy(r));
rng.simplify(2,1);
Relation AC = ConicClosure(Restrict_Range(Restrict_Domain(copy(r),copy(rng)),copy(dom)));
Relation UB = Union(copy(r),Join(copy(r),Join(AC,copy(r))));
UB.simplify(2,1);
if (detailedClosureDebug) {
fprintf(DebugFile, "UB is:\n");
UB.print_with_subs(DebugFile);
}
result = Relation::False(r);
Simple_List<Relation> firstChoice,secondChoice;
r.simplify(2,2);
Relation test = Difference(copy(r),Composition(copy(r),copy(r)));
test.simplify(2,2);
if (r.number_of_conjuncts() > test.number_of_conjuncts()) {
Relation test2 = Union(copy(test),Composition(copy(test),copy(test)));
test2.simplify(2,2);
if (Must_Be_Subset(copy(r),copy(test2))) r = test;
else if (detailedClosureDebug) {
fprintf(DebugFile, "Transitive reduction not possible:\n");
fprintf(DebugFile, "R is:\n");
r.print_with_subs(DebugFile);
fprintf(DebugFile, "test2 is:\n");
test2.print_with_subs(DebugFile);
}
}
r.make_level_carried_to(r.n_inp());
if (detailedClosureDebug) {
fprintf(DebugFile, "r is:\n");
r.print_with_subs(DebugFile);
}
for(int depth = r.n_inp(); depth >= -1; depth--)
for (DNF_Iterator d(r.query_DNF()); d.live(); d.next())
if (d.curr()->query_guaranteed_leading_0s() == depth) {
Relation C(r, d.curr());
firstChoice.append(C);
}
bool first_conj=true;
for (Simple_List_Iterator<Relation> sli(firstChoice); sli; sli++) {
if (first_conj)
first_conj=false;
else {
Relation C_plus;
bool change=TryConjunctTransitiveClosure(
copy(sli.curr()), IterationSpace, C_plus);
if (change)
sli.curr()=C_plus;
}
}
//compute closure
int maxClauses = 3+firstChoice.size()*(1+maxExpansion);
int resultConjuncts = 0;
int numFails = 0;
bool resultInexact = false;
while (!firstChoice.empty() || !secondChoice.empty()) {
Relation R_plus, R_star;
if (detailedClosureDebug) {
fprintf(DebugFile,"Main loop of TC:\n");
if (!firstChoice.empty()) {
fprintf(DebugFile,"First choice:\n");
printRelationList(firstChoice);
}
if (!secondChoice.empty()) {
fprintf(DebugFile,"Second choice:\n");
printRelationList(secondChoice);
}
}
Relation R;
if (!firstChoice.empty())
R = firstChoice.remove_front();
else R = secondChoice.remove_front();
if (detailedClosureDebug) {
fprintf(DebugFile, "Working with conjunct:\n");
R.print_with_subs(DebugFile);
}
bool known=ConjunctTransitiveClosure(copy(R),IterationSpace, R_plus, R_star);
if (!known && numFails < firstChoice.size()) {
numFails++;
firstChoice.append(R);
if (detailedClosureDebug) {
fprintf(DebugFile, "\nTry another conjunct, R is not suitable\n");
R.print_with_subs(DebugFile);
}
continue;
}
if (detailedClosureDebug) {
fprintf(DebugFile,"\nR+ is:\n");
R_plus.print_with_subs(DebugFile);
if (known) {
fprintf(DebugFile, "Known R? is :\n");
R_star.print_with_subs(DebugFile);
}
else
fprintf(DebugFile, "The R* for this relation is not calculated\n");
}
Relation R_z;
if (known) {
R_z=Difference(copy(R_star),copy(R_plus));
known = R_z.is_upper_bound_satisfiable();
if (known) {
int d = R.single_conjunct()->query_guaranteed_leading_0s();
R_z.make_level_carried_to(min(R.n_inp(),d+1));
if (R_z.query_DNF()->length() > 1) known = false;
if (detailedClosureDebug) {
fprintf(DebugFile, "\nForced R_Z to be level carried at level %d\n",min(R.n_inp(),d+1));
}
}
if (detailedClosureDebug) {
if (known) {
fprintf(DebugFile, "\nDifference between R? and R+ is:\n");
R_z.print_with_subs(DebugFile);
}
else
fprintf(DebugFile, "\nR_z is unusable\n");
}
}
else R_z = Relation::False(r);
if (!known)
numFails++;
else numFails = 0;
if (!known && numFails <= firstChoice.size()) {
firstChoice.append(R);
if (detailedClosureDebug) {
fprintf(DebugFile, "\nTry another conjunct, Rz is avaiable for R:\n");
R.print_with_subs(DebugFile);
}
continue;
}
//make N empty list
Relation N = Relation::False(r);
//append R+ to T
result = Union(result, copy(R_plus));
resultConjuncts++;
int expansion = maxClauses - (resultConjuncts + 2*firstChoice.size() + secondChoice.size());
if (expansion < 0) expansion = 0;
if (detailedClosureDebug) {
fprintf(DebugFile,"Max clauses = %d\n",maxClauses);
fprintf(DebugFile,"result conjuncts = %d\n",resultConjuncts);
fprintf(DebugFile,"firstChoice's = %d\n",firstChoice.size());
fprintf(DebugFile,"secondChoice's = %d\n",secondChoice.size());
fprintf(DebugFile,"Allowed expansion is %d\n",expansion);
}
bool firstPart=true;
if (!known && expansion == 0) {
if (detailedClosureDebug) {
fprintf(DebugFile,"Expansion = 0, R? unknown, skipping composition\n");
}
if (!resultInexact && detailedClosureDebug) fprintf(DebugFile,"RESULT BECOMES INEXACT 1\n");
resultInexact = true;
}
else
for (Simple_List_Iterator<Relation> s(firstChoice);
firstPart?
(s.live()?true:
(s = Simple_List_Iterator<Relation>(secondChoice),
firstPart = false,
s.live()))
:s.live();
s.next()) {
assert(s.live());
Relation C=(s.curr());
if (detailedClosureDebug) {
fprintf(DebugFile, "\nComposing chosen conjunct with C:\n");
C.print_with_subs(DebugFile);
}
if (!known) {
if (detailedClosureDebug) {
fprintf(DebugFile, "\nR? is unknown! No debug info here yet\n");
}
Relation C1=Composition(copy(C), copy(R_plus));
if (detailedClosureDebug) {
fprintf(DebugFile, "\nGenerating \n");
C1.print_with_subs(DebugFile);
}
C1.simplify();
Relation newStuff =
Difference(
Difference(copy(C1),copy(C)),
copy(R_plus));
newStuff.simplify();
if (detailedClosureDebug) {
fprintf(DebugFile, "New Stuff:\n");
newStuff.print_with_subs(DebugFile);
}
bool C1_contains_new_stuff = newStuff.is_upper_bound_satisfiable();
if (C1_contains_new_stuff) {
if (newStuff.has_single_conjunct())
C1 = newStuff;
if (expansion) {
N = Union(N,copy(C1));
expansion--;
}
else {
if (!resultInexact && detailedClosureDebug) fprintf(DebugFile,"RESULT BECOMES INEXACT 2\n");
resultInexact = true;
break;
}
}
else C1 = Relation::False(C1);
Relation C2(Composition(copy(R_plus),copy(C)));
if (detailedClosureDebug) {
fprintf(DebugFile, "\nGenerating \n");
C2.print_with_subs(DebugFile);
}
newStuff =
Difference(
Difference(
Difference(copy(C2),copy(C)),
copy(C1)),
copy(R_plus));
newStuff.simplify();
if (detailedClosureDebug) {
fprintf(DebugFile, "New Stuff:\n");
newStuff.print_with_subs(DebugFile);
}
if (newStuff.is_upper_bound_satisfiable()) {
if (newStuff.has_single_conjunct())
C2 = newStuff;
if (expansion) {
N = Union(N,copy(C2));
expansion--;
}
else {
if (!resultInexact && detailedClosureDebug) fprintf(DebugFile,"RESULT BECOMES INEXACT 3\n");
resultInexact = true;
break;
}
}
else C2 = Relation::False(C2);
if (C1_contains_new_stuff) {
Relation C3(Composition(copy(R_plus),copy(C1)));
if (detailedClosureDebug) {
fprintf(DebugFile, "\nGenerating \n");
C3.print_with_subs(DebugFile);
}
newStuff =
Difference(
Difference(
Difference(
Difference(copy(C3),copy(C)),
copy(C1)),
copy(C2)),
copy(R_plus));
newStuff.simplify();
if (detailedClosureDebug) {
fprintf(DebugFile, "New Stuff:\n");
newStuff.print_with_subs(DebugFile);
}
if (newStuff.is_upper_bound_satisfiable()) {
if (newStuff.has_single_conjunct())
C3 = newStuff;
if (expansion) {
N = Union(N,C3);
expansion--;
}
else {
if (!resultInexact && detailedClosureDebug) fprintf(DebugFile,"RESULT BECOMES INEXACT 4\n");
resultInexact = true;
break;
}
}
}
}
else {
Relation C_Rz(Composition(copy(C),copy(R_z)));
if (detailedClosureDebug) {
fprintf(DebugFile, "C o Rz is:\n");
C_Rz.print_with_subs(DebugFile);
}
Relation Rz_C_Rz(Composition(copy(R_z),copy(C_Rz)));
if (detailedClosureDebug) {
fprintf(DebugFile, "\nRz o C o Rz is:\n");
Rz_C_Rz.print_with_subs(DebugFile);
}
if (Equal(C,Rz_C_Rz)) {
#if defined(TC_STATS)
fprintf(statsfile,"weak test selects C?\n");
#endif
Relation tmp = Composition(C,copy(R_star));
tmp.simplify();
Relation tmp2 = Composition(copy(R_star),copy(tmp));
tmp2.simplify();
if (Must_Be_Subset(copy(tmp2),copy(tmp)))
*s = tmp;
else
*s = tmp2;
if (detailedClosureDebug) {
fprintf(DebugFile,"\nC is equal to Rz o C o Rz so R? o C o R? replaces C\n");
fprintf(DebugFile, "R? o C o R? is:\n");
(*s).print_with_subs(DebugFile);
}
}
else {
#if defined(TC_STATS)
fprintf(statsfile,"weak test fails\n");
#endif
if (Equal(C, C_Rz)) {
*s=Composition(copy(C),copy(R_star));
Relation p(Composition(copy(R_plus), copy(*s)));
p.simplify();
if (detailedClosureDebug) {
fprintf(DebugFile, "\nC is equal to C o Rz, so C o Rz replaces C\n");
fprintf (DebugFile, "C o R? is:\n");
(*s).print_with_subs(DebugFile);
fprintf (DebugFile, "R+ o C o R? is added to list N. It's :\n");
p.print_with_subs(DebugFile);
}
if (!Is_Obvious_Subset(copy(p),copy(R_plus))
&& !Is_Obvious_Subset(copy(p),copy(C))) {
if (expansion) {
p.simplify(2,2);
expansion--;
}
else {
if (!resultInexact && detailedClosureDebug) fprintf(DebugFile,"RESULT BECOMES INEXACT 5\n");
resultInexact = true;
break;
}
}
}
else {
Relation Rz_C(Composition(copy(R_z),copy(C)));
if (Equal(C,Rz_C)) {
*s=Composition(copy(R_star),copy(C));
Relation Rstar_C_Rplus(Composition(copy(*s),copy(R_plus)));
Rstar_C_Rplus.simplify();
if (detailedClosureDebug) {
fprintf(DebugFile, "\nC is equal to Rz o C , so R? o C replaces C\n");
fprintf (DebugFile, "R? o C is:\n");
(*s).print_with_subs(DebugFile);
fprintf (DebugFile, "R+ o C is added to list N. It's :\n");
Rstar_C_Rplus.print_with_subs(DebugFile);
}
if (!Is_Obvious_Subset(copy(Rstar_C_Rplus),copy(R_plus))
&& !Is_Obvious_Subset(copy(Rstar_C_Rplus),copy(C))) {
if (expansion)
N = Union(N,Rstar_C_Rplus);
else {
if (!resultInexact && detailedClosureDebug) fprintf(DebugFile,"RESULT BECOMES INEXACT 6\n");
resultInexact = true;
break;
}
}
}
else {
if (detailedClosureDebug) {
fprintf(DebugFile, "\nHave to handle it the hard way\n");
}
Relation C1=Composition(copy(C), copy(R_plus));
C1.simplify();
if (!Is_Obvious_Subset(copy(C1),copy(R_plus))
&& !Is_Obvious_Subset(copy(C1),copy(C))) {
if (expansion) {
N = Union(N,copy(C1));
expansion--;
}
else {
if (!resultInexact && detailedClosureDebug) fprintf(DebugFile,"RESULT BECOMES INEXACT 7\n");
resultInexact = true;
break;
}
}
Relation C2(Composition(copy(R_plus),copy(C)));
C2.simplify();
if (!Is_Obvious_Subset(copy(C2),copy(R_plus))
&& !Is_Obvious_Subset(copy(C2),copy(C))) {
if (expansion) {
N = Union(N,C2);
expansion--;
}
else {
if (!resultInexact && detailedClosureDebug) {
fprintf(DebugFile,"RESULT BECOMES INEXACT 8\n");
fprintf(DebugFile,"Have to discard:\n");
C2.print_with_subs(DebugFile);
}
resultInexact = true;
break;
}
}
Relation C3(Composition(copy(R_plus),C1));
C3.simplify();
if (!Is_Obvious_Subset(copy(C3),copy(R_plus)) && !Is_Obvious_Subset(copy(C3),copy(C))) {
if (expansion) {
N = Union(N,C3);
expansion--;
}
else {
if (!resultInexact && detailedClosureDebug)
fprintf(DebugFile,"RESULT BECOMES INEXACT 9\n");
resultInexact = true;
break;
}
}
}
}
}
}
}
//now we processed the first conjunct.
if (detailedClosureDebug) {
N.simplify(2,2);
fprintf(DebugFile, "\nNew conjuncts:\n");
N.print_with_subs(DebugFile);
}
N.simplify(2,2);
appendClausesToList(secondChoice,N);
}
//Did we do all conjuncts? If not, make T be inexact
result.copy_names(r);
result.simplify(2,2);
if (!result.is_exact()) {
result = Lower_Bound(result);
resultInexact = true;
}
if (resultInexact) {
Relation test(Composition(copy(result),copy(result)));
test.simplify(2,2);
if (detailedClosureDebug) {
fprintf(DebugFile, "\nResult is:\n");
result.print_with_subs(DebugFile);
fprintf(DebugFile, "\nResult composed with itself is:\n");
test.print_with_subs(DebugFile);
}
if (!Must_Be_Subset(test,copy(result))) {
result = Union(result,Intersection(UB, Relation::Unknown(result)));
}
}
#ifdef TC_STATS
{
Relation rcopy = result;
Relation test2(Composition(copy(rcopy),copy(rcopy)));
test2.simplify(2,2);
test2.remove_disjunction_with_unknown();
rcopy.remove_disjunction_with_unknown();
if (detailedClosureDebug) {
fprintf(DebugFile, "\nResult is:\n");
rcopy.print_with_subs(DebugFile);
fprintf(DebugFile, "\nResult composed with itself is:\n");
test2.print_with_subs(DebugFile);
}
if (!Must_Be_Subset(test2,copy(rcopy))) {
fprintf(statsfile,"multi TC result is inexact\n");
}
else
fprintf(statsfile,"TC result is exact%s\n", (resultInexact || !rcopy.is_exact())?" despite perceived inexactness":"");
}
#endif
#ifdef TC_STATS
}
int totalTime = clock_diff();
fprintf(statsfile, "input conj: %d output conj: %d #singe conj closures: %d time: %d\n",
in_conj, copy(result).query_DNF()->length(),
singles,
totalTime/TC_RUNS);
#endif
if (closure_presburger_debug || detailedClosureDebug) {
if (detailedClosureDebug) {
fprintf(DebugFile, "\nThe transitive closure is :\n");
result.print_with_subs(DebugFile);
}
fprintf(DebugFile, "\n\n] END Transitive closure\n\n");
}
return result;
}
Relation TransitiveClosure(NOT_CONST Relation &input_r,
int maxExpansion,
NOT_CONST Relation & input_IterationSpace) {
Relation r = consume_and_regurgitate(input_r);
Relation IterationSpace = consume_and_regurgitate(input_IterationSpace);
if (r.is_null())
return r;
if (r.n_out() == 0)
throw std::invalid_argument("transitive closure does not apply to set");
if (r.n_inp() != r.n_out())
throw std::invalid_argument("transitive closure must has the same input and output arity");
if (closure_presburger_debug) {
fprintf(DebugFile,"\nComputing Transitive closure of:\n");
r.print_with_subs(DebugFile);
fprintf(DebugFile,"\nIteration space is:\n");
IterationSpace.print_with_subs(DebugFile);
}
if (!r.is_upper_bound_satisfiable()) {
if (closure_presburger_debug)
fprintf(DebugFile, "]TC : relation is false\n");
return r;
}
Relation UB = DeltasToRelation(ConicHull(Project_Sym(Deltas(copy(r)))),
r.n_inp(),r.n_out());
if (closure_presburger_debug) {
fprintf(DebugFile,"UB is:\n");
UB.print_with_subs(DebugFile);
}
Relation conditions = Restrict_Domain(copy(UB),Domain(copy(r)));
conditions.simplify();
if (closure_presburger_debug) {
fprintf(DebugFile,"Forward reachable is:\n");
conditions.print_with_subs(DebugFile);
}
conditions = Composition(Inverse(copy(UB)),conditions);
conditions.simplify();
if (closure_presburger_debug) {
fprintf(DebugFile,"Backward/forward reachable is:\n");
conditions.print_with_subs(DebugFile);
}
conditions = Range(conditions);
conditions.simplify();
// conditions = Approximate(conditions);
// conditions.simplify();
conditions = VennDiagramForm(conditions);
conditions.simplify();
if (closure_presburger_debug) {
fprintf(DebugFile,"Condition regions are:\n");
conditions.print_with_subs(DebugFile);
}
if (conditions.is_obvious_tautology()) {
return TransitiveClosure0(r, maxExpansion, IterationSpace);
}
else {
Relation answer = Relation::False(r);
answer.copy_names(r);
answer.setup_names();
for (DNF_Iterator c(conditions.query_DNF()); c.live(); c.next()) {
Relation tmp = Relation(conditions, c.curr());
if (closure_presburger_debug) {
fprintf(DebugFile,"\nComputing Transitive closure:\n");
fprintf(DebugFile,"\nRegion:\n");
tmp.prefix_print(DebugFile);
}
Relation tmp3 = Restrict_Domain(copy(r),tmp);
tmp3.simplify(2,2);
if (closure_presburger_debug) {
fprintf(DebugFile,"\nRelation:\n");
tmp3.prefix_print(DebugFile);
}
answer = Union(answer, TransitiveClosure0(tmp3, maxExpansion,copy(IterationSpace)));
}
return answer;
}
}
/* ********************************* */
/* Function check if relation */
/* belong to d-form or */
/* uniform relaion class */
/* ********************************* */
Relation is_DForm_or_Uniform(NOT_CONST Relation &r){
Relation s = consume_and_regurgitate(r);
Relation Rtmp, Rdelta, delta;
delta = Deltas(copy(s));
Rdelta = DeltasToRelation(copy(delta), s.n_inp(), s.n_out());
Rtmp = DeltasToRelation(Project_Sym(delta), s.n_inp(), s.n_out());
Rtmp = Restrict_Domain(Rtmp, Domain(copy(Rdelta)));
Rtmp = Restrict_Range(Rtmp, Range(Rdelta));
Rdelta = copy(Rtmp);
Rtmp = Restrict_Domain(Rtmp, Domain(copy(s)));
Rtmp = Restrict_Range(Rtmp, Range(copy(s)));
if (Must_Be_Subset( copy(Rtmp), copy(s)) && \
Must_Be_Subset(copy(s), copy(Rtmp))) {
Rtmp = Relation::Null();
}
else {
Rtmp = Rdelta = Relation::Null();
}
return Rdelta;
}
/* ********************************* */
/* Get a conjunction for */
/* a given number from set */
/* of relations */
/* ********************************* */
Relation getConjunctionNr(NOT_CONST Relation &r, int conjNr) {
Relation s = consume_and_regurgitate(r);
int i = 1;
for (DNF_Iterator c(s.query_DNF()); c; c++,i++) {
if ( i == conjNr ) {
return Relation(s, c.curr());
}
}
return Relation::False(s.n_inp(), s.n_out());
}
/* ********************************* */
/* Get a common region for */
/* a given set of relations */
/* ********************************* */
Relation getCommonRegion( NOT_CONST Relation &r, const long* relTab, const long relCount) {
Relation s = consume_and_regurgitate(r);
Relation commonRegion, Rcurr;
long i = 0;
Rcurr = getConjunctionNr( copy(s), relTab[0]);
commonRegion = Union(Domain(copy(Rcurr)), Range(copy(Rcurr)));
for( i=1; i < relCount; i++ ){
Rcurr = getConjunctionNr( copy(s), relTab[i]);
commonRegion = Intersection( commonRegion, Union( Domain(copy(Rcurr)), Range(copy(Rcurr))) );
}
return commonRegion;
}
/* ********************************* */
/* Get a set of relations */
/* ********************************* */
Relation getRelationsSet( NOT_CONST Relation &r, const long* relTab, const long relCount) {
Relation s = consume_and_regurgitate(r);
Relation R = Relation::False(s.n_inp(), s.n_out());
long i = 0;
for( i=0; i < relCount; i++ ){
R = Union( R, getConjunctionNr( copy(s), relTab[i]) );
}
return R;
}
/* ********************************* */
/* Get a set of relations */
/* from a common region */
/* ********************************* */
Relation relationsOnCommonRegion( NOT_CONST Relation &r, NOT_CONST Relation ®ion ) {
Relation set = consume_and_regurgitate(r);
Relation reg = consume_and_regurgitate(region);
Relation R = Relation::True(set.n_inp(), set.n_out());
R = Restrict_Domain(R, copy(reg));
R.simplify(2,1);
R = Restrict_Range(R, reg);
R.simplify(2,1);
R = Intersection(R, set);
return R;
}
Relation compose_N(NOT_CONST Relation &input_r) {
Relation r = consume_and_regurgitate(input_r);
Relation powerR, powerR2;
r = Union(r, Identity(r.n_inp()));
powerR = copy(r);
for(;;){
if (powerR.number_of_conjuncts() > 50) {
powerR = Relation::Null();
return powerR;
}
powerR2 = Composition(copy(powerR), copy(r));
powerR2.simplify(2,1);
if (Must_Be_Subset( copy(powerR2), copy(powerR))) {
powerR2 = Relation::Null();
return powerR;
}
powerR = Relation::Null();
powerR = copy(powerR2);
powerR2 = Relation::Null();
}
}
/****************************** */
/* Check exactness of R+ */
/* */
/* Tomasz Klimek 05-06-2010 */
/****************************** */
bool checkExactness(NOT_CONST Relation &r, NOT_CONST Relation &rplus){
Relation s1 = consume_and_regurgitate(r);
Relation s2 = consume_and_regurgitate(rplus);
Relation R;
R = Composition(copy(s1), copy(s2));
R = Union(s1, R);
if( Must_Be_Subset(copy(s2), copy(R)) && \
Must_Be_Subset(copy(R), copy(s2))) {
R = Relation::Null();
s1 = Relation::Null();
return true;
}
R = Relation::Null();
s1 = Relation::Null();
return false;
}
/************************************** */
/* Calculate approximation of R* */
/* */
/* Tomasz Klimek 05-06-2010 */
/************************************** */
Relation ApproxClosure(NOT_CONST Relation &r) {
Relation s = consume_and_regurgitate(r);
Relation R = Relation::False(s.n_inp(), s.n_out());
Relation tc = Identity(s.n_inp());
Relation Rtmp;
for (DNF_Iterator c(s.query_DNF()); c; c++) {
Rtmp = Hull(Project_Sym(Deltas(Relation(s, c.curr()))), false, 1, Relation::Null());
R = Union(R, TransitiveClosure(DeltasToRelation(Rtmp,s.n_inp(),s.n_out()), 1, Relation::Null()));
}
for (DNF_Iterator c(R.query_DNF()); c; c++) {
Rtmp = Union(Identity(s.n_inp()), Relation(R, c.curr()));
tc = Composition(tc, Rtmp);
tc = Hull(tc, false, 1, Relation::Null());
}
tc = Restrict_Domain(tc,Domain(copy(s)));
tc.simplify(2,1);
tc = Restrict_Range(tc,Range(s));
tc.simplify(2,1);
tc = Intersection(tc, Relation::Unknown(tc));
return tc;
}
/************************************** */
/* Calculate R* on unbounded region */
/* */
/* Tomasz Klimek 05-06-2010 */
/************************************** */
Relation ClosureOnUnboundedRegion(NOT_CONST Relation &r) {
Relation s = consume_and_regurgitate(r);
Relation R = Relation::False(s.n_inp(), s.n_out());
Relation tc = Identity(s.n_inp());
Relation Rtmp,tcTmp;
for (DNF_Iterator c(s.query_DNF()); c; c++) {
Rtmp = is_DForm_or_Uniform(Relation(s, c.curr()));
if (!(Rtmp.is_null())) {
tcTmp = TransitiveClosure(Rtmp, 1, Relation::Null());
if (!(tcTmp.is_exact())){
tcTmp = R = Relation::Null();
/* fprintf(DebugFile,"\nTC is inexact!"); */
return tcTmp;
}
}
else {
R = Relation::Null();
/* fprintf(DebugFile,"\nR is not d-form relation!"); */
return Relation::Null();
}
R = Union(R, tcTmp);
}
for (DNF_Iterator c(R.query_DNF()); c; c++) {
Rtmp = Union(Identity(s.n_inp()), Relation(R, c.curr()));
tc = Composition(tc, Rtmp);
tc.simplify(2,1);
}
tc = Difference(tc, Identity(s.n_inp()));
return tc;
}
/******************************* */
/* Try to select sets of domain */
/* and range */
/* */
/* Tomasz Klimek 05-06-2010 */
/******************************* */
Relation SelectRegionForClosure(NOT_CONST Relation &r){
Relation s = consume_and_regurgitate(r);
Relation DR = Union(Domain(copy(s)),Range(copy(s)));
Relation region,tc,tcTmp;
region = SimpleHull(copy(DR));
region.simplify(2,1);
tc = ClosureOnUnboundedRegion(copy(s));
if (tc.is_null()) {
return tc;
}
tcTmp = Restrict_Domain(copy(tc),copy(region));
tcTmp.simplify(2,1);
tcTmp = Restrict_Range(tcTmp,region);
tcTmp.simplify(2,1);
if (checkExactness(copy(s), copy(tcTmp))) {
s = tc = Relation::Null();
return tcTmp;
}
tcTmp = Relation::Null();
region = Hull(DR,false,1,Relation::Null());
tcTmp = Restrict_Domain(copy(tc),copy(region));
tcTmp.simplify(2,1);
tcTmp = Restrict_Range(tcTmp,region);
tcTmp.simplify(2,1);
if (checkExactness(copy(s), copy(tcTmp))) {
s = tc = Relation::Null();
return tcTmp;
}
tcTmp = Relation::Null();
tc = Restrict_Domain(tc,Domain(copy(s)));
tc.simplify(2,1);
tc = Restrict_Range(tc,Domain(copy(s)));
tc.simplify(2,1);
if (checkExactness(copy(s), copy(tc))) {
s = Relation::Null();
return tc;
}
tc = Relation::Null();
return ApproxClosure(s);
}
/************************************** */
/* Calculate R* */
/* */
/* Tomasz Klimek 05-06-2010 */
/************************************** */
Relation calculateTransitiveClosure(NOT_CONST Relation &r) {
Relation s = consume_and_regurgitate(r);
Relation tc = Relation::False(s.n_inp(), s.n_out());
long* relationsSet = NULL;
Relation commonRegion, regionTmp;
Relation inputRelations;
long i,j=-1;
long N,M;
Relation R;
commonRegion = SelectRegionForClosure(copy(s));
if (commonRegion.is_null()) {
return ApproxClosure(s);
}
if (commonRegion.is_exact()) {
return commonRegion;
}
commonRegion = Relation::Null();
N = M = s.number_of_conjuncts();
relationsSet = (long*)calloc(N,sizeof(long));
if (relationsSet == NULL) {
return Relation::False(s.n_inp(), s.n_out());
}
for (; N > 1;) {
for ( i=0; i<N; i++ ) {
if ( i < j ) {
continue;
}
else if ( j == -1 ) {
relationsSet[i] = 1;
}
else if ( i > j ) {
relationsSet[i] = relationsSet[i-1] + 1;
}
else if ( i == j ) {
relationsSet[i] += 1;
}
if ( relationsSet[i] <= M ) {
j = i;
}
else {
j = i - 1;
break;
}
}
if ( j+1 == N) {
/* fprintf(DebugFile,"\n");
for(i=0;i<N;i++){
fprintf(DebugFile," %ld", relationsSet[i]);
}
fprintf(DebugFile,"\n"); */
commonRegion = getCommonRegion( copy(s), relationsSet, N);
commonRegion.simplify(2,1);
inputRelations = getRelationsSet( copy(s), relationsSet, N);
inputRelations.simplify(2,1);
/* ******************* */
/* Check on rectangle */
/* ******************* */
regionTmp = SimpleHull(copy(commonRegion));
regionTmp.simplify(2,1);
R = relationsOnCommonRegion( copy(inputRelations), regionTmp);
R.simplify(2,1);
regionTmp = SelectRegionForClosure(R);
if (regionTmp.is_exact()) {
/* fprintf(DebugFile,"\nDescribed on rectangle region\n"); */
tc = Union( tc, regionTmp );
}
else {
/* ******************* */
/* Check on hull */
/* ******************* */
R = Relation::Null();
regionTmp = Relation::Null();
regionTmp = Hull(copy(commonRegion),false,1,Relation::Null());
regionTmp.simplify(2,1);
R = relationsOnCommonRegion( copy(inputRelations), regionTmp);
R.simplify(2,1);
regionTmp = SelectRegionForClosure(R);
if (regionTmp.is_exact()) {
/* fprintf(DebugFile,"\nDescribed on Hull\n"); */
tc = Union( tc, regionTmp);
}
else {
/* ********************************** */
/* Check on sets of domain and range */
/* ********************************** */
R = Relation::Null();
regionTmp = Relation::Null();
R = relationsOnCommonRegion( copy(inputRelations), copy(commonRegion) );
R.simplify(2,1);
regionTmp = SelectRegionForClosure(R);
if (regionTmp.is_exact()) {
/* fprintf(DebugFile,"\nDescribed on sets of doamin and range\n"); */
tc = Union( tc, regionTmp );
}
else {
commonRegion = Relation::Null();
inputRelations = Relation::Null();
regionTmp = Relation::Null();
R = Relation::Null();
return ApproxClosure(s);
}
}
}
commonRegion = Relation::Null();
inputRelations = Relation::Null();
regionTmp = Relation::Null();
R = Relation::Null();
}
if ( j == -1 ) N--;
}
R = Relation::Null();
for (DNF_Iterator c(s.query_DNF()); c; c++) {
if (!Must_Be_Subset(Relation(s, c.curr()), copy(tc))) {
/* fprintf(DebugFile,"\nIs not a subset\n"); */
tc = Union( tc, SelectRegionForClosure(Relation(s, c.curr())));
}
}
if (!(tc.is_exact())){
return ApproxClosure(s);
}
tc = compose_N(tc);
if (tc.is_null()) {
return ApproxClosure(s);
}
return tc;
}
} // namespace