da: remove inconsistencies from dependence relations

[ppn.git] / eqv2.cc

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <algorithm>
#include <map>
#include <set>
#include <string>
#include <sstream>
#include <vaucanson/int_boolean_automaton.hh>

#include <isa/yaml.h>
#include <isa/pdg.h>
#include "dependence_graph.h"
#include "eqv2_options.h"
#include "msa.h"

#include "config.h"

#include <isl/constraint.h>

#ifdef HAVE_OMEGA_PLUS
#include <omega.h>

using namespace omega;

typedef std::vector<Variable_ID> varvector;
typedef std::vector<Global_Var_Decl *> globalvector;

static void create_globals(globalvector &globals, __isl_keep isl_space *dim)
{
        unsigned nparam = isl_space_dim(dim, isl_dim_param);
        for (int i = 0; i < nparam; ++i) {
                const char *name;
                name = isl_space_get_dim_name(dim, isl_dim_param, i);
                globals.push_back(new Global_Var_Decl(name));
        }
}

struct map2rel_data {
        isl_int v;

        varvector in;
        varvector out;
        varvector params;
        varvector exists;

        F_And *base;
};

static int add_constraint(__isl_take isl_constraint *c, void *user)
{
        struct map2rel_data *data = (struct map2rel_data *) user;
        Constraint_Handle h;

        if (isl_constraint_is_equality(c))
                h = data->base->add_EQ();
        else
                h = data->base->add_GEQ();
        
        for (int i = 0; i < data->in.size(); ++i) {
                isl_constraint_get_coefficient(c, isl_dim_in, i, &data->v);
                h.update_coef(data->in[i], isl_int_get_si(data->v));
        }
        for (int i = 0; i < data->out.size(); ++i) {
                isl_constraint_get_coefficient(c, isl_dim_out, i, &data->v);
                h.update_coef(data->out[i], isl_int_get_si(data->v));
        }
        for (int i = 0; i < data->params.size(); ++i) {
                isl_constraint_get_coefficient(c, isl_dim_param, i, &data->v);
                h.update_coef(data->params[i], isl_int_get_si(data->v));
        }
        for (int i = 0; i < data->exists.size(); ++i) {
                isl_constraint_get_coefficient(c, isl_dim_div, i, &data->v);
                h.update_coef(data->exists[i], isl_int_get_si(data->v));
        }
        isl_constraint_get_constant(c, &data->v);
        h.update_const(isl_int_get_si(data->v));

        isl_constraint_free(c);

        return 0;
}

static Relation basic_map2relation(isl_basic_map *bmap,
        const Variable_ID_Tuple *globals)
{
        struct map2rel_data data;

        Relation r(isl_basic_map_dim(bmap, isl_dim_in),
                   isl_basic_map_dim(bmap, isl_dim_out));
        F_Exists *e = r.add_exists();
        data.base = e->add_and();

        for (int j = 1; j <= r.n_inp(); ++j)
                data.in.push_back(r.input_var(j));
        for (int j = 1; j <= r.n_out(); ++j)
                data.out.push_back(r.output_var(j));
        for (int i = 0; i < isl_basic_map_dim(bmap, isl_dim_div); ++i)
                data.exists.push_back(e->declare());
        for (int i = 0; i < globals->size(); ++i)
                data.params.push_back(r.get_local((*globals)[i+1]));

        isl_int_init(data.v);

        isl_basic_map_foreach_constraint(bmap, &add_constraint, &data);

        isl_int_clear(data.v);

        r.finalize();
        return r;
}

static int bm2r(__isl_take isl_basic_map *bmap, void *user)
{
        Relation *r = (Relation *)user;

        *r = Union(*r, basic_map2relation(bmap, r->global_decls()));

        isl_basic_map_free(bmap);

        return 0;
}

static Relation map2relation(__isl_keep isl_map *map,
        const globalvector &globals)
{
        Relation r(isl_map_dim(map, isl_dim_in), isl_map_dim(map, isl_dim_out));
        for (int i = 0; i < globals.size(); ++i)
                r.get_local(globals[i]);
        r.add_or();
        r.finalize();

        isl_map_foreach_basic_map(map, &bm2r, &r);

        return r;
}

static void set_constraint(isl_constraint *c, Constraint_Handle h,
        varvector &in, varvector &out, varvector &params)
{
        isl_int v;

        isl_int_init(v);

        for (int i = 0; i < in.size(); ++i) {
                isl_int_set_si(v, h.get_coef(in[i]));
                isl_constraint_set_coefficient(c, isl_dim_in, i, v);
        }
        for (int i = 0; i < out.size(); ++i) {
                isl_int_set_si(v, h.get_coef(out[i]));
                isl_constraint_set_coefficient(c, isl_dim_out, i, v);
        }
        for (int i = 0; i < params.size(); ++i) {
                isl_int_set_si(v, h.get_coef(params[i]));
                isl_constraint_set_coefficient(c, isl_dim_param, i, v);
        }
        isl_int_set_si(v, h.get_const());
        isl_constraint_set_constant(c, v);

        isl_int_clear(v);
}

static void max_index(Constraint_Handle c, varvector& vv, varvector& params)
{
        for (Constr_Vars_Iter cvi(c); cvi; ++cvi) {
                Variable_ID var = (*cvi).var;
                if (find(vv.begin(), vv.end(), var) == vv.end() &&
                    find(params.begin(), params.end(), var) == params.end())
                        vv.push_back(var);
        }
}

static void set_constraint(isl_constraint *c, Constraint_Handle h,
                                varvector &vv, varvector &params)
{
        isl_int v;

        isl_int_init(v);

        for (int i = 0; i < vv.size(); ++i) {
                isl_int_set_si(v, h.get_coef(vv[i]));
                isl_constraint_set_coefficient(c, isl_dim_set, i, v);
        }
        for (int i = 0; i < params.size(); ++i) {
                isl_int_set_si(v, h.get_coef(params[i]));
                isl_constraint_set_coefficient(c, isl_dim_param, i, v);
        }
        isl_int_set_si(v, h.get_const());
        isl_constraint_set_constant(c, v);

        isl_int_clear(v);
}

static isl_basic_set *disjunct2basic_set(DNF_Iterator &di, isl_space *dim,
                                            varvector &vv, varvector &params)
{
        int d = vv.size();
        isl_constraint *c;
        isl_basic_set *bset;
        isl_local_space *ls;

        for (EQ_Iterator ei = (*di)->EQs(); ei; ++ei)
                max_index((*ei), vv, params);
        for (GEQ_Iterator gi = (*di)->GEQs(); gi; ++gi)
                max_index((*gi), vv, params);

        dim = isl_space_add_dims(dim, isl_dim_set, vv.size() - d);
        bset = isl_basic_set_universe(isl_space_copy(dim));
        ls = isl_local_space_from_space(dim);

        for (EQ_Iterator ei = (*di)->EQs(); ei; ++ei) {
                c = isl_equality_alloc(isl_local_space_copy(ls));
                set_constraint(c, (*ei), vv, params);
                bset = isl_basic_set_add_constraint(bset, c);
        }
        for (GEQ_Iterator gi = (*di)->GEQs(); gi; ++gi) {
                c = isl_inequality_alloc(isl_local_space_copy(ls));
                set_constraint(c, (*gi), vv, params);
                bset = isl_basic_set_add_constraint(bset, c);
        }

        isl_local_space_free(ls);

        bset = isl_basic_set_project_out(bset, isl_dim_set, d, vv.size() - d);
        vv.resize(d);

        return bset;
}

static __isl_give isl_map *relation2map(__isl_take isl_space *dim, Relation &r)
{
        varvector vv;
        varvector params;
        isl_ctx *ctx;
        isl_space *dim_set;
        isl_map *map;

        for (int j = 1; j <= r.n_inp(); ++j)
                vv.push_back(r.input_var(j));
        for (int j = 1; j <= r.n_out(); ++j)
                vv.push_back(r.output_var(j));

        const Variable_ID_Tuple *globals = r.global_decls();
        for (int i = 0; i < globals->size(); ++i)
                params.push_back(r.get_local((*globals)[i+1]));

        ctx = isl_space_get_ctx(dim);
        dim_set = isl_space_set_alloc(ctx, params.size(), vv.size());
        map = isl_map_empty(isl_space_copy(dim));

        for (DNF_Iterator di(r.query_DNF()); di; ++di) {
                isl_basic_set *bset;
                isl_basic_map *bmap;

                bset = disjunct2basic_set(di, isl_space_copy(dim_set), vv, params);
                bmap = isl_basic_map_from_basic_set(bset, isl_space_copy(dim));
                map = isl_map_union(map, isl_map_from_basic_map(bmap));
        }

        isl_space_free(dim_set);
        isl_space_free(dim);

        return map;
}

/* Compute transitive closure using Omega+.
 * Note that the result computed by Omega also contains the identity
 * mapping on the intersection of domain and range, but we will include
 * the entire identity mapping later, so this shouldn't be a problem.
 */
static __isl_give isl_map *omega_closure(__isl_take isl_map *map)
{
        isl_ctx *ctx = isl_map_get_ctx(map);
        struct options *options = isl_ctx_peek_eqv2_options(ctx);
        isl_space *dim;
        globalvector globals;

        dim = isl_map_get_space(map);
        create_globals(globals, dim);

        Relation R = map2relation(map, globals);
        isl_map_free(map);

        Relation TC;
        if (options->try_lower) {
                int exact;
                TC = TransitiveClosure(copy(R));
                Relation test = Composition(copy(TC), copy(R));

                exact = TC.is_exact() &&
                        Must_Be_Subset(copy(R), copy(TC)) &&
                        Must_Be_Subset(test, copy(TC));
                if (!exact)
                        TC = Upper_Bound(ApproxClosure(R));
        } else
                TC = Upper_Bound(ApproxClosure(R));

        map = relation2map(dim, TC);

        return map;
}
#endif

/* Compute transitive closure */
static __isl_give isl_map *closure(__isl_take isl_map *map)
{
        isl_ctx *ctx = isl_map_get_ctx(map);
        struct options *options = isl_ctx_peek_eqv2_options(ctx);

#ifdef HAVE_OMEGA_PLUS
        if (options->omega)
                return omega_closure(map);
#endif

        return isl_map_transitive_closure(map, NULL);
}

using namespace vcsn::int_boolean_automaton;
using namespace vcsn;

template <class Exp_, class Dispatch_>
class acc_visitor : public algebra::KRatExpMatcher<
      acc_visitor<Exp_, Dispatch_>, Exp_, isl_map *, Dispatch_ >
{
public:
        typedef typename Exp_::monoid_elt_value_t       monoid_elt_value_t;
        typedef acc_visitor<Exp_, Dispatch_>            this_class;
        typedef isl_map *                               return_type;

private:
        const MSA &msa;

public:
        acc_visitor(const MSA &msa) : msa(msa) {}

        INHERIT_CONSTRUCTORS(this_class, Exp_, isl_map *, Dispatch_);

        MATCH__(Product, lhs, rhs) {
                return isl_map_apply_range(match(lhs), match(rhs));
        } END

        MATCH__(Sum, lhs, rhs) {
                return isl_map_union(match(lhs), match(rhs));
        } END

        MATCH_(Star, node) {
                isl_space *dim;
                isl_map *map;
                map = closure(match(node));
                dim = isl_map_get_space(map);
                return isl_map_union(map, isl_map_identity(dim));
        } END

        MATCH__(LeftWeight, w, node) {
                assert(0);
                return 0;
        } END

        MATCH__(RightWeight, node, w) {
                assert(0);
                return 0;
        } END

        MATCH_(Constant, m) {
                isl_map *map;
                typename monoid_elt_value_t::const_iterator i = m.begin();
                assert(i != m.end());
                map = isl_map_copy(msa.firing[*i]);
                ++i;
                assert(i == m.end());

                return map;
        } END

        /* By construction, there should always be a path */
        MATCH(Zero) {
                assert(0);
                return 0;
        } END

        /* There can be no zero-length path */
        MATCH(One) {
                assert(0);
                return 0;
        } END
};

template<typename Exp>
isl_map *accessibility(const MSA &msa, const Exp &kexp)
{
        acc_visitor<Exp, algebra::DispatchFunction<Exp> > m(msa);
        return m.match(kexp);
}

static unsigned update_out_dim(pdg::PDG *pdg, unsigned out_dim)
{
        for (int i = 0; i < pdg->arrays.size(); ++i) {
                pdg::array *array = pdg->arrays[i];
                if (array->type != pdg::array::output)
                        continue;
                if (array->dims.size() + 1 > out_dim)
                        out_dim = array->dims.size() + 1;
        }

        return out_dim;
}

/* The input arrays of the two programs are supposed to be the same,
 * so they should at least have the same dimension.  Make sure
 * this is true, because we depend on it later on.
 */
static int check_input_arrays(dependence_graph *dg1, dependence_graph *dg2)
{
        for (int i = 0; i < dg1->input_computations.size(); ++i)
                for (int j = 0; j < dg2->input_computations.size(); ++j) {
                        if (strcmp(dg1->input_computations[i]->operation,
                                   dg2->input_computations[j]->operation))
                                continue;
                        if (dg1->input_computations[i]->dim ==
                            dg2->input_computations[j]->dim)
                                continue;
                        fprintf(stderr,
                            "input arrays \"%s\" do not have the same dimension\n",
                            dg1->input_computations[i]->operation);
                        return -1;
                }

        return 0;
}

/* Compute a map between the domains of the output nodes of
 * the dependence graphs that expresses that array1 in dg1
 * and array2 in dg2 are equal.
 * This map is then returned as a set.
 */
static __isl_give isl_set *compute_equal_array(
        dependence_graph *dg1, dependence_graph *dg2, int array1, int array2)
{
        const char *array_name = dg1->output_arrays[array1];
        isl_set *out1 = isl_set_copy(dg1->out->domain);
        isl_set *out2 = isl_set_copy(dg2->out->domain);
        unsigned dim1 = isl_set_n_dim(out1);
        unsigned dim2 = isl_set_n_dim(out2);

        out1 = isl_set_fix_dim_si(out1, dim1 - 1, array1);
        out2 = isl_set_fix_dim_si(out2, dim2 - 1, array2);
        out1 = isl_set_remove_dims(out1, isl_dim_set, dim1 - 1, 1);
        out2 = isl_set_remove_dims(out2, isl_dim_set, dim2 - 1, 1);
        int equal = isl_set_is_equal(out1, out2);
        isl_set_free(out1);
        isl_set_free(out2);
        assert(equal >= 0);
        if (!equal) {
                fprintf(stderr, "different output domains for array %s\n",
                                array_name);
                return NULL;
        }

        out1 = isl_set_copy(dg1->out->domain);
        out2 = isl_set_copy(dg2->out->domain);
        out1 = isl_set_fix_dim_si(out1, dim1 - 1, array1);
        out2 = isl_set_fix_dim_si(out2, dim2 - 1, array2);
        out1 = isl_set_coalesce(out1);
        out2 = isl_set_coalesce(out2);
        isl_space *dim = isl_space_map_from_set(isl_set_get_space(out2));
        isl_map *id = isl_map_identity(dim);
        id = isl_map_remove_dims(id, isl_dim_in, isl_map_n_in(id) - 1, 1);
        id = isl_map_apply_domain(id, isl_map_copy(id));
        id = isl_map_intersect_domain(id, out1);
        id = isl_map_intersect_range(id, out2);

        return isl_map_wrap(id);
}

/* Compute a map between the domains of the output nodes of
 * the dependence graphs that expresses that arrays with
 * the same names in both graphs are equal.
 * The map is returned as a set.
 */
static __isl_give isl_set *compute_equal_output(
        dependence_graph *dg1, dependence_graph *dg2)
{
        int i1 = 0, i2 = 0;
        isl_space *dim;
        isl_map *map;
        isl_set *equal_output;
        unsigned d;

        dim = isl_set_get_space(dg1->out->domain);
        d = isl_space_dim(dim, isl_dim_set);
        map = isl_map_empty(isl_space_map_from_set(dim));
        equal_output = isl_map_wrap(map);

        while (i1 < dg1->output_arrays.size() || i2 < dg2->output_arrays.size()) {
                int cmp;
                cmp = i1 == dg1->output_arrays.size() ? 1 :
                      i2 == dg2->output_arrays.size() ? -1 :
                      strcmp(dg1->output_arrays[i1], dg2->output_arrays[i2]);
                if (cmp < 0) {
                        fprintf(stdout,
                                "Array \"%s\" only in program 1; not checked\n",
                                dg1->output_arrays[i1]);
                        ++i1;
                } else if (cmp > 0) {
                        fprintf(stdout,
                                "Array \"%s\" only in program 2; not checked\n",
                                dg2->output_arrays[i2]);
                        ++i2;
                } else {
                        isl_set *eq_array;
                        eq_array = compute_equal_array(dg1, dg2, i1, i2);
                        equal_output = isl_set_union(equal_output, eq_array);
                        ++i1;
                        ++i2;
                }
        }

        equal_output = isl_set_set_tuple_name(equal_output, "s0");

        return equal_output;
}

static int check_equivalence(isl_ctx *ctx,
        dependence_graph *dg1, dependence_graph *dg2)
{
        int root;
        MSA msa;
        isl_printer *p;
        unsigned d;
        int ok = 1;

        if (check_input_arrays(dg1, dg2))
                return -1;

        p = isl_printer_to_file(ctx, stdout);

        root = msa.add(dg1->out, dg2->out);

        alphabet_t alpha;
        for (int i = 0; i < msa.firing.size(); ++i)
                alpha.insert(i);

        automaton_t a = make_automaton(alpha);
        std::vector<hstate_t> state;
        for (int i = 0; i < msa.state.size(); ++i)
                state.push_back(a.add_state());

        for (int i = 0; i < msa.transition.size(); ++i) {
                int s1 = msa.transition[i].first.first;
                int s2 = msa.transition[i].first.second;
                int f = msa.transition[i].second;
                a.add_letter_transition(state[s1], state[s2], f);
        }

        a.set_initial(state[root]);

        isl_set *equal_output;
        equal_output = compute_equal_output(dg1, dg2);

        for (int i = 0; i < msa.input_state.size(); ++i) {
                isl_map *acc;
                isl_set *set;
                isl_map *map;
                isl_map *id;

                a.set_final(state[msa.input_state[i]]);
                rat_exp_t ratexp = aut_to_exp(a);
                a.clear_final();

                acc = accessibility(msa, ratexp.value());
                set = isl_set_apply(isl_set_copy(equal_output),
                                        isl_map_copy(acc));
                map = isl_set_unwrap(set);

                id = isl_map_identity(isl_map_get_space(map));
                map = isl_map_subtract(map, id);

                set = isl_map_wrap(map);
                set = isl_set_set_tuple_name(set,
                                    isl_map_get_tuple_name(acc, isl_dim_out));
                set = isl_set_apply(set, isl_map_reverse(acc));
                d = isl_set_dim(set, isl_dim_set) / 2;
                set = isl_set_remove_dims(set, isl_dim_set, d, d);
                if (!isl_set_is_empty(set)) {
                        ok = 0;
                        isl_printer_print_str(p, "NOT equivalent");
                        isl_printer_end_line(p);
                        set = isl_set_coalesce(set);
                        isl_printer_print_set(p, set);
                        isl_printer_end_line(p);
                }

                isl_set_free(set);
        }

        /* Consider each fail state in turn as the corresponding
         * domains may have different dimensions.
         */
        for (int i = 0; i < msa.fail_state.size(); ++i) {
                isl_map *acc;
                isl_set *set;

                a.set_final(state[msa.fail_state[i]]);
                rat_exp_t ratexp = aut_to_exp(a);
                a.clear_final();

                acc = accessibility(msa, ratexp.value());
                acc = isl_map_intersect_domain(acc, isl_set_copy(equal_output));

                set = isl_map_domain(acc);
                d = isl_set_dim(set, isl_dim_set) / 2;
                set = isl_set_remove_dims(set, isl_dim_set, d, d);
                if (!isl_set_is_empty(set)) {
                        ok = 0;
                        isl_printer_print_str(p, "NOT equivalent");
                        isl_printer_end_line(p);
                        set = isl_set_coalesce(set);
                        isl_printer_print_set(p, set);
                        isl_printer_end_line(p);
                }

                isl_set_free(set);
        }

        isl_set_free(equal_output);

        if (!ok)
                isl_printer_print_str(p, "NOT ");
        isl_printer_print_str(p, "OK");
        isl_printer_end_line(p);

        isl_printer_free(p);

        return 0;
}

int main(int argc, char *argv[])
{
        struct isl_ctx *ctx;
        struct options *options = eqv2_options_new_with_defaults();
        isl_set *context = NULL;
        dependence_graph *dg[2];

        argc = eqv2_options_parse(options, argc, argv, ISL_ARG_ALL);

        ctx = isl_ctx_alloc_with_options(&options_args, options);
        if (!ctx) {
                fprintf(stderr, "Unable to allocate ctx\n");
                return -1;
        }

        if (options->context)
                context = isl_set_read_from_str(ctx, options->context);

        pdg::PDG *pdg[2];
        unsigned out_dim = 0;
        for (int i = 0; i < 2; ++i) {
                FILE *in;
                in = fopen(options->program[i], "r");
                if (!in) {
                        fprintf(stderr, "Unable to open %s\n", options->program[i]);
                        return -1;
                }
                pdg[i] = yaml::Load<pdg::PDG>(in, ctx);
                fclose(in);
                if (!pdg[i]) {
                        fprintf(stderr, "Unable to read %s\n", options->program[i]);
                        return -1;
                }
                out_dim = update_out_dim(pdg[i], out_dim);

                if (context &&
                    pdg[i]->params.size() != isl_set_dim(context, isl_dim_param)) {
                        fprintf(stdout,
                            "Parameter dimension mismatch; context ignored\n");
                        isl_set_free(context);
                        context = NULL;
                }
        }

        for (int i = 0; i < 2; ++i) {
                dg[i] = pdg_to_dg(pdg[i], out_dim, isl_set_copy(context));
                pdg[i]->free();
                delete pdg[i];
        }

        int res = check_equivalence(ctx, dg[0], dg[1]);

        for (int i = 0; i < 2; ++i)
                delete dg[i];
        isl_set_free(context);
        isl_ctx_free(ctx);

        return res;
}