#include <omega.h>
#include <omega/Relations.h>
#include <basic/DynamicArray.h>
#include <omega/reach.h>

namespace omega {

typedef DynamicArray1<Relation> Rel_Array1;
typedef DynamicArray2<Relation> Rel_Array2;

// This is from parallelism.c, modified

static void closure_rel(Rel_Array2 &trans, int i, int k, int j) {
  Relation tik;

  if (trans[k][k].is_upper_bound_satisfiable()) {
    Relation tkk = TransitiveClosure(copy(trans[k][k]));
    tkk.simplify(2,4);
    tik = Composition(tkk, copy(trans[i][k]));
    tik.simplify(2,4);
    tik = Union(copy(trans[i][k]), tik);
    tik.simplify(2,4);
  }
  else {
    tik = trans[i][k];
  }
  Relation fresh;
  Relation tkj = trans[k][j];
  fresh = Composition(tkj, tik);
  fresh.simplify(2,4);
  trans[i][j] = Union(trans[i][j], fresh);
  trans[i][j].simplify(2,4);

#if 0
  fprintf(DebugFile, "%d -> %d -> %d\n", i, k, j);
  trans[i][j].print_with_subs(DebugFile);
#endif
}


static void close_rels(Rel_Array2 &trans,int n_nodes) {
  for (int k=1; k<=n_nodes; k++)
    for (int i=1; i<=n_nodes; i++)
      if (trans[i][k].is_upper_bound_satisfiable())
        for (int j=1; j<=n_nodes; j++)
          if (trans[k][j].is_upper_bound_satisfiable())
            closure_rel(trans, i, k, j);
}


void dump_rels(Rel_Array2 &a, reachable_information *reachable_info) {
  int i,j;
  int n_nodes = reachable_info->node_names.size();
  for(i = 1; i <= n_nodes; i++)
    for(j = 1; j <= n_nodes; j++) {
      fprintf(stderr,"t[%s][%s] = ",
              (reachable_info->node_names[i]).c_str(),
              (reachable_info->node_names[j]).c_str());
      a[i][j].print_with_subs(stderr);
    }
}


void dump_sets(Rel_Array1 &a, reachable_information *reachable_info) {
  int i;
  int n_nodes = reachable_info->node_names.size();
  for(i = 1; i <= n_nodes; i++) {
    fprintf(stderr,"r[%s] = ", (reachable_info->node_names[i]).c_str());
    a[i].print_with_subs(stderr);
  }
}



Rel_Array1 *Reachable_Nodes(reachable_information *reachable_info) {
  Tuple<std::string> &node_names = reachable_info->node_names;
  Tuple<int> &arity = reachable_info->node_arity;
  Rel_Array2 &transitions = reachable_info->transitions;
  Rel_Array1 &start_nodes = reachable_info->start_nodes;

  int n_nodes = node_names.size(),i,j;

#define DUMP_INITIAL 1
#define DUMP_CLOSED 1

  if(DUMP_INITIAL && relation_debug){
    fprintf(stderr,"Initially:\n");
    dump_rels(transitions, reachable_info);
  }

  close_rels(transitions,n_nodes);

  if(DUMP_CLOSED && relation_debug) {
    fprintf(stderr,"Closed:\n");
    dump_rels(transitions, reachable_info);
  }

  Rel_Array1 *finalp = 
    new Rel_Array1("node");
  Rel_Array1 &final = *finalp;
  final.resize(n_nodes+1);
  for (i=1; i<=n_nodes; i++)
    final[i] = Relation::False(arity[i]);

  for(i = 1; i <= n_nodes; i++)
    for(j = 1; j <= n_nodes; j++)
      if(start_nodes[i].is_upper_bound_satisfiable())
        final[j] = Union(final[j],
                         Composition(copy(transitions[i][j]),
                                     copy(start_nodes[i])));
  return finalp; 
}

static void compute_initially_reachable(Rel_Array1 &r,
                                        Rel_Array1 &start_nodes,
                                        Rel_Array2 &,
                                        Rel_Array2 &closed,
                                        Rel_Array1 &end_nodes,
                                        int n_nodes, Tuple<int> &arity){
  for(int n = 1; n <= n_nodes; n++)
    r[n] = Relation::False(arity[n]);

  for(int i = 1; i <= n_nodes; i++)
    for(int j = 1; j <= n_nodes; j++)
      r[i] = Union(r[i],
                   Range(Restrict_Domain(
                           Restrict_Range(copy(closed[j][i]),
                                          copy(end_nodes[i])),
                           copy(start_nodes[j]))));
}


static bool iterate(Rel_Array1 &r, Rel_Array2 &, Rel_Array2 &closed,
                    Rel_Array1 &, int n_nodes) {
  bool changed;

  changed = false;
  for(int j = 1; j <= n_nodes; j++) {
    for(int i = 1; i <= n_nodes; i++) {
      /* look for additional steps from interesting states */
      Relation new_rj =  Range(Restrict_Domain(copy(closed[i][j]),
                                               copy(r[i])));
      if(!Must_Be_Subset(copy(new_rj),copy(r[j]))) {
        r[j] = Union(r[j],new_rj);
        r[j].simplify(2,2);
        changed = true;
      }
    }
  }
  return changed;
}




Rel_Array1 *I_Reachable_Nodes(reachable_information *reachable_info) {
  bool changed;

  Tuple<std::string> &node_names = reachable_info->node_names;
  int n_nodes = node_names.size();
  Tuple<int> &arity = reachable_info->node_arity;
  Rel_Array2 &transitions = reachable_info->transitions;
  Rel_Array1 &start_nodes = reachable_info->start_nodes;
  Rel_Array2 closed("node number","node number");
  closed.resize(n_nodes+1,n_nodes+1);  // abuse of dynamic arrays
  Rel_Array1 *rp = new Rel_Array1("node number");
  Rel_Array1 &r = *rp;
  r.resize(n_nodes+1);  // abuse of dynamic arrays

  int i,j;

  Rel_Array1 end_nodes("Hi!");
  end_nodes.resize(n_nodes+1); // for future use
  for(int n = 1; n <= n_nodes; n++) end_nodes[n] = Relation::True(arity[n]);

  for(j = 1; j <= n_nodes; j++) {
    closed[j][j] = TransitiveClosure(copy(transitions[j][j]));
    for(i = 1; i <= n_nodes; i++)
      if (i != j)
        closed[i][j] = transitions[i][j];
  }

  compute_initially_reachable(r,start_nodes,transitions,closed,end_nodes,
                              n_nodes,arity);

#define DUMP_INITIAL 1
#define DUMP_CLOSED 1

  if(DUMP_INITIAL && relation_debug > 1) {
    fprintf(stderr,"Closed:\n");
    dump_rels(closed, reachable_info);
  }
  if(DUMP_INITIAL && relation_debug) {
    fprintf(stderr,"start nodes:\n");
    dump_sets(start_nodes, reachable_info);
    fprintf(stderr,"Initially reachable:\n");
    dump_sets(r, reachable_info);
  }

  changed = true;
  int iterations = 0, max_iterations = 1000;
  while(changed && iterations < max_iterations) {
    changed = iterate(r,transitions,closed,start_nodes,n_nodes);
    iterations++;
  }
  if(relation_debug)
    fprintf(stdout,"[Iterations to convergence: %d]\n",iterations);
  return rp;
}

} // namespace