diff --git a/.clang-tidy b/.clang-tidy
index 794cf1c2c2e5776c2fe499a05bf0f70f74281967..ff6e7f6449a967ba36e50b247349c8b34371ad50 100644
--- a/.clang-tidy
+++ b/.clang-tidy
@@ -1,3 +1,4 @@
 Checks: 'cppcoreguidelines-*,google-*,modernize-use-nullptr,performance-*,portability-*,readability-*,-cppcoreguidelines-owning-memory'
-WarningsAsErrors: true
+WarningsAsErrors: 'true'
 FormatStyle: google
+MinimumVariableNameLength: 1
diff --git a/README.md b/README.md
index 84c2fea761824130c78ede7a7fb3deb78a23e5ee..5ab184bba29d6f52cbe511a66c01098a9412594c 100644
--- a/README.md
+++ b/README.md
@@ -7,6 +7,8 @@ generates the abstract paths.
 
 - `gcc >= 9.3.0`.
 - `cmake`
+- `flex`
+- `bison`
 
 ## Build
 
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 9919fda2a9f83c1ae0ea0061fdf400a5c0db3d9f..4f395b1d24f14320b086cb99ef312652b8fbb496 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -3,9 +3,9 @@ cmake_minimum_required(VERSION 3.5 FATAL_ERROR)
 
 # main program
 add_executable(sogMBT main.cpp
-  MDGraph.cpp
-  RdPBDD.cpp
-  Modular_Obs_Graph.cpp
+        sog.cpp
+        petri_net_sog.cpp
+        modular_sog.cpp
 )
 
 # link libraries
diff --git a/src/Class_of_state.hpp b/src/Class_of_state.hpp
deleted file mode 100644
index 21c51231bbae577f5375501b89dd60aaac080ea0..0000000000000000000000000000000000000000
--- a/src/Class_of_state.hpp
+++ /dev/null
@@ -1,26 +0,0 @@
-#ifndef CLASS_OF_STATE
-#define CLASS_OF_STATE
-using namespace std;
-#include "bdd.h"
-// #include <map.h>
-#include <set>
-#include <unordered_map>
-#include <vector>
-// #include <pair>
-// #include <ext/hash_map>
-typedef set<int> Set;
-class Class_Of_State {
- public:
-  Class_Of_State() { boucle = blocage = Visited = 0; }
-  Set firable;
-  bdd class_state;
-  bool boucle;
-  bool blocage;
-  bool Visited;
-  void* Class_Appartenance;
-  vector<pair<Class_Of_State*, int> > Predecessors, Successors;
-  pair<Class_Of_State*, int> LastEdge;
-};
-typedef pair<Class_Of_State*, int> Edge;
-typedef vector<Edge> Edges;
-#endif
diff --git a/src/MDGraph.cpp b/src/MDGraph.cpp
deleted file mode 100644
index eacf70c9afb861c0ac3d82f2014b6adc52f20985..0000000000000000000000000000000000000000
--- a/src/MDGraph.cpp
+++ /dev/null
@@ -1,158 +0,0 @@
-
-/********              Graph.cpp     *******/
-#include "MDGraph.hpp"
-
-#include <fstream>
-
-bdd *Tab;
-/*********                  setInitialState    *****/
-
-void MDGraph::setInitialState(Class_Of_State *c) {
-  currentstate = initialstate = c;
-}
-/*----------------------find()----------------*/
-Class_Of_State *MDGraph::find(Class_Of_State *c) {
-  for (MetaGrapheNodes::const_iterator i = GONodes.begin();
-       !(i == GONodes.end()); i++)
-    // if((c->class_state.id()==(*i)->class_state.id())&&(c->blocage==(*i)->blocage)&&(c->boucle==(*i)->boucle))
-    if (c->class_state.id() == (*i)->class_state.id()) return *i;
-  return NULL;
-}
-
-/*----------------------insert() ------------*/
-void MDGraph::insert(Class_Of_State *c) {
-  c->Visited = false;
-  this->GONodes.push_back(c);
-  nbStates++;
-  // cout<<"on ajoute ag => nb agregats= "<<nbStates<<endl;
-}
-
-/*----------------------NbBddNod()------------------------*/
-int MDGraph::NbBddNode(Class_Of_State *S, size_t &nb) {
-  if (S->Visited == false) {
-    // cout<<"insertion du meta etat numero :"<<nb<<"son id est
-    // :"<<S->class_state.id()<<endl; cout<<"sa taille est
-    // :"<<bdd_nodecount(S->class_state)<<" noeuds \n";
-    Tab[nb - 1] = S->class_state;
-    S->Visited = true;
-    int bddnode = bdd_nodecount(S->class_state);
-    int size_succ = 0;
-    for (Edges::const_iterator i = S->Successors.begin();
-         !(i == S->Successors.end()); i++) {
-      if ((*i).first->Visited == false) {
-        nb++;
-        size_succ += NbBddNode((*i).first, nb);
-      }
-    }
-    return size_succ + bddnode;
-
-  } else
-    return 0;
-}
-
-/*----------------------Visualisation du graphe------------------------*/
-void MDGraph::printCompleteInformation() {
-  cout << "\n\nGRAPH SIZE : \n";
-  cout << "\n\tNB MARKING : " << nbMarking;
-  cout << "\n\tNB NODES : " << nbStates;
-  cout << "\n\tNB ARCS : " << nbArcs << endl;
-  // cout<<" \n\nCOMPLETE INFORMATION ?(y/n)\n";
-  // char c;
-  // cin>>c;
-  // InitVisit(initialstate,n);
-  Tab = new bdd[(int)nbStates];
-  size_t n = 1;
-  // cout<<"NB BDD NODE : "<<NbBddNode(initialstate,n)<<endl;
-  NbBddNode(initialstate, n);
-  cout << "NB BDD NODE : " << bdd_anodecount(Tab, (int)nbStates) << endl;
-  // cout<<"Shared Nodes : "<<bdd_anodecount(Tab,nbStates)<<endl;
-  InitVisit(initialstate, 1);
-  // int toto;cin>>toto;
-  // bdd Union=UnionMetaState(initialstate,1);
-  // cout<<"a titre indicatif taille de l'union : "<<bdd_nodecount(Union)<<endl;
-  /*if(c=='y'||c=='Y')
-  {
-          size_t n=1;
-           printGraph(initialstate,n);
-  }*/
-}
-/*----------------------InitVisit()------------------------*/
-void MDGraph::InitVisit(Class_Of_State *S, size_t nb) {
-  if (nb <= nbStates) {
-    S->Visited = false;
-    for (Edges::const_iterator i = S->Successors.begin();
-         !(i == S->Successors.end()); i++) {
-      if ((*i).first->Visited == true) {
-        nb++;
-        InitVisit((*i).first, nb);
-      }
-    }
-  }
-}
-/*********                  printGraph    *****/
-
-void MDGraph::printGraph(Class_Of_State *s, size_t &nb) {
-  if (nb <= nbStates) {
-    cout << "\nSTATE NUMBER " << nb << " : \n";
-    s->Visited = true;
-    printsuccessors(s);
-    getchar();
-    printpredecessors(s);
-    getchar();
-    Edges::const_iterator i;
-    for (i = s->Successors.begin(); !(i == s->Successors.end()); i++) {
-      if ((*i).first->Visited == false) {
-        nb++;
-        printGraph((*i).first, nb);
-      }
-    }
-  }
-}
-
-/*----------------------generate the file of the reachability graph
- * -------------------*/
-void MDGraph::generate_Dotfile_of_reachability_graph(string filenamep) {
-  string filename("./result/" + filenamep + ".dot");
-  fstream file;
-  cout << "le graph contient : " << GONodes.size() << endl;
-  file.open(filename, std::ios_base::out);
-  file << "digraph reachab_graph {" << endl;
-  size_t n = 1;
-  for (auto i : GONodes) {
-    for (auto k : i->Successors) {
-      file << "ag_" << i->class_state.id() << " -> "
-           << "ag_" << k.first->class_state.id() << " [ label =  \"t"
-           << k.second + 1 << "\"]" << endl;
-    }
-  }
-
-  file << "ag_" << initialstate->class_state.id() << "[initialstate]" << endl;
-  file << " }" << endl;
-}
-
-/*---------void print_successors_class(Class_Of_State *)------------*/
-void MDGraph::printsuccessors(Class_Of_State *s) {
-  Edges::const_iterator i;
-  cout << bddtable << s->class_state << endl;
-  if (s->boucle) cout << "\n\tON BOUCLE DESSUS AVEC EPSILON\n";
-  if (s->blocage) cout << "\n\tEXISTENCE D'UN ETAT BLOCANT\n";
-  cout << "\n\tSES SUCCESSEURS SONT  ( " << s->Successors.size() << " ) :\n\n";
-  getchar();
-  for (i = s->Successors.begin(); !(i == s->Successors.end()); i++) {
-    cout << " \t- t" << (*i).second << " ->";
-    cout << bddtable << (*i).first->class_state << endl;
-    getchar();
-  }
-}
-/*---------void printpredescessors(Class_Of_State *)------------*/
-void MDGraph::printpredecessors(Class_Of_State *s) {
-  Edges::const_iterator i;
-  cout << "\n\tSES PREDESCESSEURS SONT  ( " << s->Predecessors.size()
-       << " ) :\n\n";
-  getchar();
-  for (i = s->Predecessors.begin(); !(i == s->Predecessors.end()); i++) {
-    cout << " \t- t" << (*i).second << " ->";
-    cout << bddtable << (*i).first->class_state << endl;
-    getchar();
-  }
-}
diff --git a/src/MDGraph.hpp b/src/MDGraph.hpp
deleted file mode 100644
index 7be028bf66a138080d6cb2c9561a79c40994d218..0000000000000000000000000000000000000000
--- a/src/MDGraph.hpp
+++ /dev/null
@@ -1,45 +0,0 @@
-
-
-/******************    Graph.hpp  **************************/
-
-#ifndef _MDGRAPH_
-
-#define _MDGRAPH_
-#include <iostream>
-#include <list>
-#include <string>
-#include <vector>
-
-using namespace std;
-
-#include "Class_of_state.hpp"
-typedef set<int> Set;
-typedef vector<Class_Of_State *> MetaGrapheNodes;
-class MDGraph {
- private:
-  void printGraph(Class_Of_State *, size_t &);
-
- public:
-  MetaGrapheNodes GONodes;
-  void Reset();
-  Class_Of_State *initialstate;
-  Class_Of_State *currentstate;
-  double nbStates;
-  double nbMarking;
-  double nbArcs;
-  Class_Of_State *find(Class_Of_State *);
-  Edges &get_successor(Class_Of_State *);
-  void printsuccessors(Class_Of_State *);
-  int NbBddNode(Class_Of_State *, size_t &);
-  void InitVisit(Class_Of_State *S, size_t nb);
-  void printpredecessors(Class_Of_State *);
-  void addArc() { nbArcs++; }
-  void insert(Class_Of_State *);
-  MDGraph() { nbStates = nbArcs = nbMarking = 0; }
-  void setInitialState(
-      Class_Of_State *);  // Define the initial state of this graph
-  void printCompleteInformation();
-  void generate_Dotfile_of_reachability_graph(
-      string filenamep);  // as input for Graphwalker
-};
-#endif
diff --git a/src/Modular_Class_of_state.hpp b/src/Modular_Class_of_state.hpp
deleted file mode 100644
index 25591e1a07b31fd8239fab77e715d801806ab603..0000000000000000000000000000000000000000
--- a/src/Modular_Class_of_state.hpp
+++ /dev/null
@@ -1,57 +0,0 @@
-#ifndef MODULAR_CLASS_OF_STATE
-#define MODULAR_CLASS_OF_STATE
-// #include<hash_map.h>
-#include <set>
-
-#include "Net.hpp"
-// #include <vector>
-#include "bdd.h"
-//_____________________   Vecteur d'entier : repr�sente les id des bdd relatifs
-// � un �tat synchronis� de n GO _________________________ typedef vector<bdd>
-// SynchronizedState; struct less_SynchronizedState
-//{
-// bool operator()(const SynchronizedState,const SynchronizedState)const;
-//};
-// inline bool less_SynchronizedState::operator ()(const SynchronizedState
-// p1,const SynchronizedState p2)const
-//{
-
-// for(int i=0;i<p1.size();i++)
-//	if(p1[i].id()>p2[i].id())
-//		return false;
-// return true;
-
-//}
-// typedef set<SynchronizedState, less_SynchronizedState> Gen_Prod_Synch;
-class Modular_Class_Of_State {
-  // friend struct less_modular_class_of_state;
- public:
-  Modular_Class_Of_State() { boucle = blocage = Visited = 0; }
-  vector<bdd> State;
-  bool boucle;
-  bool blocage;
-  bool Visited;
-  set<pair<Modular_Class_Of_State *, string> > Predecessors;
-  set<pair<Modular_Class_Of_State *, string> > Successors;
-  friend ostream &operator<<(ostream &os, const Modular_Class_Of_State &c) {
-    os << "{";
-    for (unsigned int k = 0; k < c.State.size(); k++)
-      os << c.State[k].id() << ", ";
-    os << "}";
-    return (os);
-  }
-};
-/*____________________ MODULAR EDGE________________*/
-class Modular_Class_Of_State;
-typedef pair<Modular_Class_Of_State *, string> Modular_Edge;
-struct LessModularEdge {
-  bool operator()(const Modular_Edge, const Modular_Edge) const;
-};
-inline bool LessModularEdge::operator()(const Modular_Edge a1,
-                                        const Modular_Edge a2) const {
-  // cout<<"on compare les etiquettes des arcs : "<<a1.second<<" avec
-  // "<<a2.second<<endl;
-  return (a1.first < a2.first);
-}
-typedef set<Modular_Edge, LessModularEdge> Modular_Edges;
-#endif
diff --git a/src/Modular_Obs_Graph.cpp b/src/Modular_Obs_Graph.cpp
deleted file mode 100644
index caa80dcc10beb87bb03ac2b9db39fc5dff5132c4..0000000000000000000000000000000000000000
--- a/src/Modular_Obs_Graph.cpp
+++ /dev/null
@@ -1,181 +0,0 @@
-
-/********              Graph.cpp     *******/
-#include "Modular_Obs_Graph.hpp"
-// #include<conio.h>
-bdd *Temp;  // tableau interm�diaire pour calculer la taille (nb bdd) du graphe
-/*********                  setInitialState    *****/
-
-void Modular_Obs_Graph::setInitialState(Modular_Class_Of_State *c) {
-  currentstate = initialstate = c;
-}
-/*----------------------find()----------------*/
-Modular_Class_Of_State *Modular_Obs_Graph::find(Modular_Class_Of_State *c) {
-  bool arret;
-  for (Modular_Obs_Graph_Nodes::const_iterator i = GONodes.begin();
-       !(i == GONodes.end()); i++) {
-    // if((c->blocage!=(*i)->blocage)||(c->boucle!=(*i)->boucle))
-    //{
-    arret = false;
-    for (unsigned int k = 0; ((k < (c->State).size()) && (!arret)); k++)
-      if (!(c->State[k] == (*i)->State[k])) arret = true;
-    if (!arret) return *i;
-    //}
-  }
-  return NULL;
-}
-/*----------------------find2()----------------*/
-// version modulaire on v�rifie la projection
-Modular_Class_Of_State *Modular_Obs_Graph::find2(Modular_Class_Of_State *c,
-                                                 Set SRConcernes) {
-  bool arret;
-  for (Modular_Obs_Graph_Nodes::const_iterator i = GONodes.begin();
-       !(i == GONodes.end()); i++) {
-    // if((c->blocage!=(*i)->blocage)||(c->boucle!=(*i)->boucle))
-    //{
-    arret = false;
-    for (unsigned int k = 0; ((k < (c->State).size()) && (!arret)); k++)
-      if (!(SRConcernes.find(k) == SRConcernes.end()))
-        if (!(c->State[k] == (*i)->State[k])) arret = true;
-    if (!arret && (c->blocage == (*i)->blocage) && (c->boucle == (*i)->boucle))
-      return *i;
-    //}
-  }
-  return NULL;
-}
-/*-----------------------AddARc2----------------*/
-void Modular_Obs_Graph::addArc(Modular_Class_Of_State *Pred,
-                               Modular_Class_Of_State *Suc, const char *t) {
-  // cout<<"ici addArc entre "<<*Pred<<" et "<<*Suc<<endl;
-  Modular_Edge arc = Modular_Edge(Suc, t);
-  Modular_Edge cra = Modular_Edge(Pred, t);
-  if (Pred->Successors.find(arc) == Pred->Successors.end()) {
-    // cout<<"OK \n";
-    Pred->Successors.insert(Pred->Successors.begin(), arc);
-    Suc->Predecessors.insert(Suc->Predecessors.begin(), cra);
-    nbArcs++;
-  }
-  // else
-  //	cout<<"KO \n";
-  // int toto;cin>>toto;
-}
-/*----------------------insert() ------------*/
-void Modular_Obs_Graph::insert(Modular_Class_Of_State *c) {
-  c->Visited = false;
-  this->GONodes.push_back(c);
-  nbStates++;
-}
-
-/*----------------------InitVisit()------------------------*/
-void Modular_Obs_Graph::InitVisit(Modular_Class_Of_State *S, size_t nb) {
-  // cout<<"____________ nb = __________________ "<<nb<<endl;
-  // cout<<"____________ nbStates = __________________ "<<nbStates<<endl;
-  if (nb <= nbStates) {
-    //	cout<<"nb < = nbStates\n";
-    S->Visited = false;
-    for (Modular_Edges::const_iterator i = S->Successors.begin();
-         !(i == S->Successors.end()); i++) {
-      if ((*i).first->Visited == true) {
-        nb++;
-        InitVisit((*i).first, nb);
-      }
-    }
-  }
-}
-/*----------------------NbBddNod()------------------------*/
-void Modular_Obs_Graph::TAB_BDDNODES(Modular_Class_Of_State *S, size_t &nb) {
-  if (S->Visited == false) {
-    // cout<<"ETAT NON VISITE : ";
-    // cout<<*S<<endl;
-    // int tt;cin>>tt;
-    for (unsigned int k = 0; k < S->State.size(); k++, nb++)
-      Temp[nb - 1] = S->State[k];
-    nb--;
-    S->Visited = true;
-    // int bddnode=0;
-    // for(k=0;k<S->State.size();k++)
-    //	bddnode+=bdd_nodecount(S->State[k]);
-    // int size_succ=0;
-    for (Modular_Edges::const_iterator i = S->Successors.begin();
-         !(i == S->Successors.end()); i++) {
-      if ((*i).first->Visited == false) {
-        nb++;
-        TAB_BDDNODES((*i).first, nb);
-      }
-    }
-  }
-}
-/*------------------------------------------Affichage du graphe -------------*/
-void Modular_Obs_Graph::printCompleteInformation(int nbsubnets) {
-  cout << "\n\nGRAPH SIZE : \n";
-  cout << "\n\tNB MARKING : " << nbMarking;
-  cout << "\n\tNB NODES : " << nbStates;
-  cout << "\n\tNB ARCS : " << nbArcs << endl;
-  cout << " \n\nCOMPLETE INFORMATION ?(y/n)\n";
-  char c;
-  cin >> c;
-  // InitVisit(initialstate,n);
-  Temp = new bdd[nbStates * nbsubnets];
-  size_t n = 1;
-  // cout<<"NB BDD NODE : "<<NbBddNode(initialstate,n)<<endl;
-  TAB_BDDNODES(initialstate, n);
-  cout << "NB BDD NODE : " << bdd_anodecount(Temp, nbStates * nbsubnets)
-       << endl;
-  // cout<<"Shared Nodes : "<<bdd_anodecount(Tab,nbStates)<<endl;
-  InitVisit(initialstate, 1);
-  if (c == 'y' || c == 'Y') {
-    size_t n = 1;
-    printGraph(initialstate, n);
-  }
-}
-/*********                  printGraph    *****/
-
-void Modular_Obs_Graph::printGraph(Modular_Class_Of_State *s, size_t &nb) {
-  // cout<<"Print Graph \n";
-  // int toto;cin>>toto;
-  if (nb <= nbStates) {
-    cout << "\nSTATE NUMBER " << nb << " : \n";
-    s->Visited = true;
-    printsuccessors(s);
-    getchar();
-    printpredecessors(s);
-    getchar();
-    Modular_Edges::const_iterator i;
-    for (i = s->Successors.begin(); !(i == s->Successors.end()); i++) {
-      if ((*i).first->Visited == false) {
-        nb++;
-        printGraph((*i).first, nb);
-      }
-    }
-  }
-}
-/*-------------------------Reset()----------------------------*/
-void Modular_Obs_Graph::Reset() {
-  currentstate = NULL;
-  nbArcs = nbMarking = nbStates = 0;
-}
-/*---------void print_successors_class(Class_Of_State *)------------*/
-void Modular_Obs_Graph::printsuccessors(Modular_Class_Of_State *s) {
-  Modular_Edges::const_iterator i;
-  cout << *s << endl;
-  if (s->boucle) cout << "\n\tON BOUCLE DESSUS AVEC EPSILON\n";
-  if (s->blocage) cout << "\n\tEXISTENCE D'UN ETAT BLOCANT\n";
-  cout << "\n\tSES SUCCESSEURS SONT  ( " << s->Successors.size() << " ) :\n\n";
-  getchar();
-  for (i = s->Successors.begin(); !(i == s->Successors.end()); i++) {
-    cout << " \t---" << (*i).second << " -->";
-    cout << *((*i).first) << endl;
-    getchar();
-  }
-}
-/*---------void printpredescessors(Class_Of_State *)------------*/
-void Modular_Obs_Graph::printpredecessors(Modular_Class_Of_State *s) {
-  Modular_Edges::const_iterator i;
-  cout << "\n\tSES PREDESCESSEURS SONT  ( " << s->Predecessors.size()
-       << " ) :\n\n";
-  getchar();
-  for (i = s->Predecessors.begin(); !(i == s->Predecessors.end()); i++) {
-    cout << " \t---" << (*i).second << " -->";
-    cout << *(*i).first << endl;
-    getchar();
-  }
-}
diff --git a/src/Modular_Obs_Graph.hpp b/src/Modular_Obs_Graph.hpp
deleted file mode 100644
index ea11ae32b9b0aa97dd56d5da764b1ba5096363cd..0000000000000000000000000000000000000000
--- a/src/Modular_Obs_Graph.hpp
+++ /dev/null
@@ -1,44 +0,0 @@
-
-
-/******************    Graph.hpp  **************************/
-
-#ifndef _Modular_ObsGraph_
-
-#define _Modular_ObsGraph_
-#include <time.h>
-
-#include <iostream>
-#include <vector>
-
-#include "Modular_Class_of_state.hpp"
-typedef set<int> Set;
-typedef vector<Modular_Class_Of_State *> Modular_Obs_Graph_Nodes;
-class Modular_Obs_Graph {
- private:
-  void printGraph(Modular_Class_Of_State *, size_t &);
-  Modular_Obs_Graph_Nodes GONodes;
-
- public:
-  void Reset();
-  Modular_Class_Of_State *initialstate;
-  Modular_Class_Of_State *currentstate;
-  size_t nbStates;
-  size_t nbMarking;
-  size_t nbArcs;
-  void Liberer(Modular_Class_Of_State *S);
-  Modular_Class_Of_State *find(Modular_Class_Of_State *);
-  Modular_Class_Of_State *find2(Modular_Class_Of_State *, Set);
-  Modular_Edges &get_successor(Modular_Class_Of_State *);
-  void printsuccessors(Modular_Class_Of_State *);
-  void InitVisit(Modular_Class_Of_State *S, size_t nb);
-  void TAB_BDDNODES(Modular_Class_Of_State *, size_t &);
-  void printpredecessors(Modular_Class_Of_State *);
-  void addArc() { nbArcs++; }
-  void addArc(Modular_Class_Of_State *, Modular_Class_Of_State *, const char *);
-  void insert(Modular_Class_Of_State *);
-  Modular_Obs_Graph() { nbStates = nbArcs = nbMarking = 0; }
-  void setInitialState(
-      Modular_Class_Of_State *);  // Define the initial state of this graph
-  void printCompleteInformation(int nbsubnets);
-};
-#endif
diff --git a/src/RdPBDD.cpp b/src/RdPBDD.cpp
deleted file mode 100644
index 1d6ada5e18e072d75a5a9ca144b1a0973d78fee1..0000000000000000000000000000000000000000
--- a/src/RdPBDD.cpp
+++ /dev/null
@@ -1,1153 +0,0 @@
-/* -*- C++ -*- */
-#include <iostream>
-#include <map>
-#include <set>
-#include <stack>
-#include <string>
-#include <unordered_map>
-#include <vector>
-
-#include "Net.hpp"
-using namespace std;
-#include <math.h>
-
-#include "RdPBDD.hpp"
-#include "bvec.h"
-int NbIt;
-int itext, itint;
-int MaxIntBdd;
-bdd *TabMeta;
-int nbmetastate;
-double old_size;
-const vector<class Place> *vplaces = NULL;
-void my_error_handler(int errcode) {
-  // cout<<"errcode = "<<errcode<<endl;
-  if (errcode == BDD_RANGE) {
-    // Value out of range : increase the size of the variables...
-    // but which one???
-    bdd_default_errhandler(errcode);
-  } else {
-    bdd_default_errhandler(errcode);
-  }
-}
-
-/*****************************************************************/
-/*                         printhandler                          */
-/*****************************************************************/
-void printhandler(ostream &o, int var) {
-  o << (*vplaces)[var / 2].name;
-  if (var % 2) o << "_p";
-}
-
-/*****************************************************************/
-/*                         class Trans                           */
-/*****************************************************************/
-Trans::Trans(const bdd &v, bddPair *p, const bdd &r, const bdd &pr,
-             const bdd &pre, const bdd &post)
-    : var(v), pair(p), rel(r), prerel(pr), Precond(pre), Postcond(post) {}
-// Franchissement avant
-bdd Trans::operator()(const bdd &n) const {
-  bdd res = bdd_relprod(n, rel, var);
-  res = bdd_replace(res, pair);
-  return res;
-}
-// Franchissement arri�re
-bdd Trans::operator[](const bdd &n) const {
-  bdd res = bdd_relprod(n, prerel, var);
-  res = bdd_replace(res, pair);
-  return res;
-}
-
-/*****************************************************************/
-/*                         class RdPBDD                          */
-/*****************************************************************/
-
-RdPBDD::RdPBDD(const net &R, map<int, int> observables, Set NonObservables,
-               int BOUND, bool init) {
-  int nbPlaces = R.places.size(), i, domain;
-  int nbTransitions = R.transitions.size();
-  vector<Place>::const_iterator p;
-
-  bvec *v = new bvec[nbPlaces];
-  bvec *vp = new bvec[nbPlaces];
-  bvec *prec = new bvec[nbPlaces];
-  bvec *postc = new bvec[nbPlaces];
-  int *idv = new int[nbPlaces];
-  int *idvp = new int[nbPlaces];
-  int *nbbit = new int[nbPlaces];
-  if (!init) bdd_init(1000000, 1000000);
-  // the error handler
-  bdd_error_hook((bddinthandler)my_error_handler);
-  //_______________
-  transitions = R.transitions;
-  for (auto i : observables) {
-    Observable.insert(i.first);
-  };
-
-  NonObservable = NonObservables;
-  transitionName = R.transitionName;
-  Nb_places = R.places.size();
-  // cout<<"Nombre de places : "<<Nb_places<<endl;
-  // cout<<"Nombre de transitions : "<<nbTransitions<<endl;
-  // cout<<"Derniere place : "<<R.places[Nb_places-1].name<<endl;
-  /* place domain, place bvect, place initial marking and place name */
-  // domains
-  i = 0;
-  for (i = 0, p = R.places.begin(); p != R.places.end(); i++, p++) {
-    if (p->hasCapacity()) {
-      domain = p->capacity + 1;
-    } else {
-      domain = BOUND + 1;  // the default domain
-    }
-    // variables are created one by one (implying contigue binary variables)
-    fdd_extdomain(&domain, 1);
-    // cout<<"nb created var : "<<bdd_varnum()<<endl;
-    fdd_extdomain(&domain, 1);
-    // cout<<"nb created var : "<<bdd_varnum()<<endl;
-  }
-
-  // bvec
-  currentvar = bdd_true();
-  for (i = 0; i < nbPlaces; i++) {
-    nbbit[i] = fdd_varnum(2 * i);
-    // cout<<"nb var pour "<<2*i<<" = "<<fdd_domainsize(2*i)<<endl;
-    v[i] = bvec_varfdd(2 * i);
-    vp[i] = bvec_varfdd(2 * i + 1);
-    // cout<<"nb var pour "<<2*i+1<<" = "<<fdd_domainsize(2*i+1)<<endl;
-    currentvar = currentvar & fdd_ithset(2 * i);
-  }
-
-  // initial marking
-  M0 = bdd_true();
-  for (i = 0, p = R.places.begin(); p != R.places.end(); i++, p++)
-    M0 = M0 & fdd_ithvar(2 * i, p->marking);
-
-  // place names
-  vplaces = &R.places;
-  fdd_strm_hook(printhandler);
-
-  /* Transition relation */
-  for (vector<Transition>::const_iterator t = R.transitions.begin();
-       t != R.transitions.end(); t++) {
-    int np = 0;
-    bdd rel = bdd_true(), var = bdd_true(), prerel = bdd_true();
-    bdd Precond = bdd_true(), Postcond = bdd_true();
-    bddPair *p = bdd_newpair();
-
-    for (i = 0; i < nbPlaces; i++) {
-      prec[i] = bvec_con(nbbit[i], 0);
-      postc[i] = bvec_con(nbbit[i], 0);
-    }
-    set<int> adjacentPlace;
-
-    // arcs pre
-    for (vector<pair<int, int>>::const_iterator it = t->pre.begin();
-         it != t->pre.end(); it++) {
-      adjacentPlace.insert(it->first);
-      prec[it->first] =
-          prec[it->first] + bvec_con(nbbit[it->first], it->second);
-    }
-    // arcs post
-    for (vector<pair<int, int>>::const_iterator it = t->post.begin();
-         it != t->post.end(); it++) {
-      adjacentPlace.insert(it->first);
-      postc[it->first] =
-          postc[it->first] + bvec_con(nbbit[it->first], it->second);
-    }
-
-    np = 0;
-    for (set<int>::const_iterator it = adjacentPlace.begin();
-         it != adjacentPlace.end(); it++) {
-      idv[np] = 2 * (*it);
-      idvp[np] = 2 * (*it) + 1;
-      var = var & fdd_ithset(2 * (*it));
-      // Image
-      //  precondition
-      rel = rel & (v[*it] >= prec[*it]);
-      Precond = Precond & (v[*it] >= prec[*it]);
-      // postcondition
-      rel = rel & (vp[*it] == (v[*it] - prec[*it] + postc[*it]));
-      // Pre image __________
-      // precondition
-      prerel = prerel & (v[*it] >= postc[*it]);
-      // postcondition
-      Postcond = Postcond & (v[*it] >= postc[*it]);
-      prerel = prerel & (vp[*it] == (v[*it] - postc[*it] + prec[*it]));
-      //___________________
-      // capacit�
-      if (R.places[*it].hasCapacity())
-        rel = rel & (vp[*it] <= bvec_con(nbbit[*it], R.places[*it].capacity));
-      np++;
-    }
-    fdd_setpairs(p, idvp, idv, np);
-    relation.push_back(Trans(var, p, rel, prerel, Precond, Postcond));
-  }
-  delete[] v;
-  delete[] vp;
-  delete[] prec;
-  delete[] postc;
-  delete[] idv;
-  delete[] idvp;
-  delete[] nbbit;
-}
-
-/*******************************************************************************************/
-/*----------------------------------- Reachability space using bdd ----------*/
-bdd RdPBDD::ReachableBDD1() {
-  bdd M1;
-  bdd M2 = M0;
-  double d, tps;
-  d = (double)clock() / (double)CLOCKS_PER_SEC;
-  NbIt = 0;
-  MaxIntBdd = bdd_nodecount(M0);
-  while (M1 != M2) {
-    M1 = M2;
-    for (vector<Trans>::const_iterator i = relation.begin();
-         i != relation.end(); i++) {
-      bdd succ = (*i)(M2);
-      M2 = succ | M2;
-      // cout << bdd_nodecount(succ) << endl;
-      // if(succ!=bddfalse)
-    }
-    NbIt++;
-    int Current_size = bdd_nodecount(M2);
-    if (MaxIntBdd < Current_size) MaxIntBdd = Current_size;
-    // cout<<"Iteration numero "<<NbIt<<" nb node de reached
-    // :"<<bdd_nodecount(M2)<<endl;
-    //		//cout << bdd_nodecount(M2) << endl;
-  }
-  // cout << endl;
-  tps = ((double)(clock()) / CLOCKS_PER_SEC) - d;
-  // cout<<"-----------------------------------------------------\n";
-  // cout << "CONSTRUCTION TIME  " << tps << endl;
-  // cout<<"Max Intermediary BDD "<<MaxIntBdd<<endl;
-  // cout<<"Nombre d'iteration : "<<NbIt<<endl;
-  // bdd Cycle=EmersonLey(M2,0);
-  // cout<<Cycle.id()<<endl;
-  return M2;
-}
-bdd RdPBDD::ReachableBDD2() {
-  bdd M1;
-  bdd M2 = M0;
-  double d, tps;
-  d = (double)clock() / (double)CLOCKS_PER_SEC;
-  NbIt = 0;
-  MaxIntBdd = bdd_nodecount(M0);
-  while (M1 != M2) {
-    M1 = M2;
-    bdd Reached;
-    for (vector<Trans>::const_iterator i = relation.begin();
-         i != relation.end(); i++) {
-      bdd succ = (*i)(M2);
-      Reached = succ | Reached;
-      // cout << bdd_nodecount(succ) << endl;
-      // if(succ!=bddfalse)
-    }
-    NbIt++;
-    M2 = M2 | Reached;
-    int Current_size = bdd_nodecount(M2);
-    if (MaxIntBdd < Current_size) MaxIntBdd = Current_size;
-    // cout<<"Iteration numero "<<NbIt<<" nb node de reached
-    // :"<<bdd_nodecount(M2)<<endl;
-    //		//cout << bdd_nodecount(M2) << endl;
-  }
-  // cout << endl;
-  tps = ((double)(clock()) / CLOCKS_PER_SEC) - d;
-  // cout<<"-----------------------------------------------------\n";
-  // cout << "CONSTRUCTION TIME  " << tps << endl;
-  // cout<<"Max Intermediary BDD "<<MaxIntBdd<<endl;
-  // cout<<"Nombre d'iteration : "<<NbIt<<endl;
-  return M2;
-}
-bdd RdPBDD::ReachableBDD3() {
-  double d, tps;
-  d = (double)clock() / (double)CLOCKS_PER_SEC;
-  bdd New, Reached, From;
-  Reached = From = M0;
-  NbIt = 0;
-  do {
-    NbIt++;
-    bdd succ;
-    for (vector<Trans>::const_iterator i = relation.begin();
-         i != relation.end(); i++)
-      succ = (*i)(From) | succ;
-    New = succ - Reached;
-    From = New;
-    Reached = Reached | New;
-    // cout<<"Iteration numero "<<NbIt<<" nb node de reached
-    // :"<<bdd_nodecount(Reached)<<endl;
-  } while (New != bddfalse);
-  tps = (double)clock() / (double)CLOCKS_PER_SEC - d;
-  // cout << "TPS CONSTRUCTION : "<<tps<<endl;
-  return Reached;
-}
-/*----------------Fermeture transitive sur les transitions non observ�es ---*/
-bdd RdPBDD::Accessible_epsilon2(bdd Init) {
-  bdd Reached, New, From;
-  Reached = From = Init;
-  do {
-    bdd succ;
-    for (Set::const_iterator i = NonObservable.begin();
-         !(i == NonObservable.end()); i++)
-
-      succ = relation[(*i)](From) | succ;
-    New = succ - Reached;
-    From = New;
-    Reached = Reached | New;
-  } while (New != bddfalse);
-  cout << endl;
-  return Reached;
-}
-bdd RdPBDD::Accessible_epsilon(bdd From) {
-  bdd M1;
-  bdd M2 = From;
-  // int it=0;
-  // cout<<"Ici accessible epsilon"<<endl;
-  do {
-    M1 = M2;
-    for (Set::const_iterator i = NonObservable.begin();
-         !(i == NonObservable.end()); i++) {
-      // cout<<" unobs : "<<transitions[*i].name<<endl;
-      bdd succ = relation[(*i)](M2);
-      M2 = succ | M2;
-    }
-    // TabMeta[nbmetastate]=M2;
-    // int intsize=bdd_anodecount(TabMeta,nbmetastate+1);
-    // if(MaxIntBdd<intsize)
-    //    MaxIntBdd=intsize;
-    // it++;
-    //	cout << bdd_nodecount(M2) << endl;
-  } while (M1 != M2);
-  // cout << endl;
-  //  cout<<"Fin accessible epsilon"<<endl;
-  return M2;
-} /*------------------------  StepForward()  --------------*/
-bdd RdPBDD::StepForward2(bdd From) {
-  // //cout<<"Debut Step Forward \n";
-  bdd Res;
-  for (Set::const_iterator i = NonObservable.begin();
-       !(i == NonObservable.end()); i++) {
-    bdd succ = relation[(*i)](From);
-    Res = Res | succ;
-  }
-  // cout<<"Fin Step Forward \n";
-  return Res;
-}
-bdd RdPBDD::StepForward(bdd From) {
-  // cout<<"Debut Step Forward \n";
-  bdd Res = From;
-  for (Set::const_iterator i = NonObservable.begin();
-       !(i == NonObservable.end()); i++) {
-    bdd succ = relation[(*i)](Res);
-    Res = Res | succ;
-  }
-  // cout<<"Fin Step Forward \n";
-  return Res;
-}
-/*------------------------  StepBackward2()  --------------*/
-bdd RdPBDD::StepBackward2(bdd From) {
-  bdd Res;
-  // cout<<"Ici Step Back : From.id() = "<<From.id()<<endl;
-  for (vector<Trans>::const_iterator t = relation.begin(); t != relation.end();
-       t++) {
-    bdd succ = (*t)[From];
-    Res = Res | succ;
-  }
-  // cout<<"Res.id() = "<<Res.id()<<endl;
-  return Res;
-}
-/*------------------------  StepBackward()  --------------*/
-bdd RdPBDD::StepBackward(bdd From) {
-  bdd Res = From;
-  for (vector<Trans>::const_iterator t = relation.begin(); t != relation.end();
-       t++) {
-    bdd succ = (*t)[Res];
-    Res = Res | succ;
-  }
-  return Res;
-}
-
-pair<int, bdd> RdPBDD::StepBackward1(bdd From, Class_Of_State *agr) {
-  pair<int, bdd> p;
-  for (auto t : NonObservable) {
-    // cout<<"on traite transition "<<t+1<<endl;
-    bdd succ = (relation[t])[From];
-    // cout<<"on aura le bdd  "<<succ<<endl;
-
-    // la fonction qui retourne le bdd precedent avec la transition t
-
-    if ((succ != bdd_false()) & ((succ &= agr->class_state) != bdd_false())) {
-      // cout<<"non null et appartient a agr !  "<<endl;
-      p.first = t;
-      p.second = succ;
-      break;
-    }
-  }
-
-  return p;
-}
-
-bool operator<(const pair<bdd, bdd> r, const pair<bdd, bdd> s) {
-  return r.first.id() == s.first.id();
-}
-
-/*---------------------------GetSuccessor()-----------*/
-bdd RdPBDD::get_successor(bdd From, int t) { return relation[t](From); }
-/*-------------------------Firable Obs()--------------*/
-Set RdPBDD::firable_obs(bdd State) {
-  Set res;
-  for (Set::const_iterator i = Observable.begin(); !(i == Observable.end());
-       i++) {
-    {
-      bdd succ = relation[(*i)](State);
-      if (succ != bddfalse) res.insert(*i);
-    }
-  }
-  return res;
-}
-
-/**-----------------------les points de sortie d'un agr avec une transition
- * donn�e t ---------*/
-
-void RdPBDD::bdd_firable_obs(Class_Of_State *agr, int t) {
-  /*Set res;
-  for(Set::const_iterator i=Observable.begin();!(i==Observable.end());i++)
-  {
-        {bdd succ =  relation[(*i)](State);
-        if(succ!=bddfalse)
-          res.insert(*i);}
-   }
-  return res;*/
-  // cout<<"les differents bdd de l agregat "<<agr->class_state<<endl;
-}
-
-/*---------------------------------------------------------*/
-
-int findMax(set<int> my_set) {
-  // Get the maximum element
-  int max_element;
-  if (!my_set.empty()) max_element = *(my_set.rbegin());
-
-  // return the maximum element
-  return max_element;
-}
-
-int choisir(set<int> cov, set<int> fire) {
-  std::set<int>::iterator it = fire.begin();
-
-  while (it != fire.end()) {
-    if (const bool is_in = cov.find(*it) == cov.end()) {
-      return *it;
-    }
-    it++;
-  }
-
-  return -1;
-}
-
-void reinit_cycle(vector<int> sw, map<int, int> &tranobs) {
-  for (auto i : sw) {
-    if (tranobs[i] > 0) {
-      tranobs[i] = 2;
-    }
-  }
-}
-
-/*---------------------------------extraction des chemins
- * ---------------------------------*/
-typedef vector<int> chem;
-
-set<chem> RdPBDD::chem_obs(MDGraph &g, map<int, int> tranobs) {
-  set<chem> cto;
-  vector<int> sw;
-  bool ancien = true;
-  set<int> cov;
-  typedef pair<Class_Of_State *, bdd> couple;
-  typedef pair<couple, Set> Pair;
-  typedef stack<Pair> pile;
-  pile st;
-  int t;
-
-  double d, tps;
-  d = (double)clock() / (double)CLOCKS_PER_SEC;
-  TabMeta = new bdd[100000];
-  nbmetastate = 0;
-  MaxIntBdd = 0;
-  NbIt = 0;
-  itext = itint = 0;
-  Class_Of_State *reached_class;
-  Set fire;
-
-  Class_Of_State *c = new Class_Of_State;  // construction duu premier agregat
-  {
-    bdd Complete_meta_state = Accessible_epsilon(M0);
-    fire = firable_obs(Complete_meta_state);
-    c->class_state = Complete_meta_state;
-    TabMeta[nbmetastate] = c->class_state;
-    nbmetastate++;
-    old_size = bdd_nodecount(c->class_state);
-    st.push(Pair(couple(c, Complete_meta_state), fire));
-  }
-  g.setInitialState(c);  // init
-  g.insert(c);
-  // //cout << "premier agregat construit " << endl;
-  g.nbMarking += bdd_pathcount(c->class_state);
-  /* for (auto f: fire) {
-      //cout << "le nb de fires de"<<c->class_state.id()<<" est " << fire.size()
-   << " : t" << f+1 << endl;
-   }*/
-  while (!st.empty()) {
-    // //cout<<"il reste des elmts dans st " <<endl;
-    Pair e = st.top();
-    st.pop();
-    nbmetastate--;
-
-    if (!e.second.empty())  // firable non vide
-    {
-      // //cout<<"il reste des fire" <<endl;
-      // choisir transition des transitions firables
-      int t = choisir(cov, e.second);
-      if (t != -1) {
-        ancien = false;
-        tranobs[t]--;
-        sw.push_back(t);
-        if (tranobs[t] == 0) {
-          cov.insert(t);
-          if (cov.size() == Observable.size()) {
-            if (!ancien) {
-              cto.insert(sw);
-            }
-            return cto;
-          }
-        }
-      }
-
-      else {
-        t = *e.second.begin();
-        sw.push_back(t);
-      }
-      // //cout<<"+++++++++++++++on a trouv� la transition t"<<t+1<<endl;
-      e.second.erase(t);
-
-      // //cout << "on ajoute t"<<t+1 <<" au chemin"<< endl;
-
-      // //cout << "ici on couvre  la transition t" << t+1 <<"qui a un second de
-      // "<<tranobs[t]<< " et le nb de trans couv = " << cov.size()
-      // <<"/"<<Observable.size()<< endl;
-      st.push(e);
-
-      double nbnode;
-      reached_class = new Class_Of_State;
-      {
-        // //cout << "on traite un bdd avac la trans t" <<t+1<< endl;
-        bdd Complete_meta_state =
-            Accessible_epsilon(get_successor(e.first.second, t));
-        // cout << "on crre un bdd" << endl;
-        reached_class->class_state = Complete_meta_state;
-        Class_Of_State *pos = g.find(reached_class);
-        nbnode = bdd_pathcount(reached_class->class_state);
-
-        if (!pos)  // nouvel agregat
-        {
-          // //cout << "nouvel agregat" << endl;
-          fire = firable_obs(Complete_meta_state);
-          /* for (auto f: fire) {
-             //cout << "le nb de fires de  "<<reached_class->class_state.id()<<"
-           est " << fire.size() << " : t" << f+1 << endl;
-           }*/
-          // st.push(Pair(couple(reached_class, Complete_meta_state), fire));
-          nbmetastate++;
-          e.first.first->Successors.insert(e.first.first->Successors.begin(),
-                                           Edge(reached_class, t));
-          reached_class->Predecessors.insert(
-              reached_class->Predecessors.begin(), Edge(e.first.first, t));
-          // //cout<<"les pred sont avec t"<<Edge(e.first.first,
-          // t).second+1<<endl;
-          g.addArc();
-          g.insert(reached_class);  // ajouter l'agregat
-          if (fire.empty()) {
-            // cout << "pas de successeurs " << endl;
-
-            if (!ancien) {
-              cto.insert(sw);
-            }
-            reinit_cycle(sw, tranobs);
-
-            sw.pop_back();
-            // cout << "on enleve du chemin"<< endl;
-
-          } else {
-            // cout << "il y a  de successeurs " << endl;
-            st.push(Pair(couple(reached_class, Complete_meta_state), fire));
-          }
-
-        } else {  // agregat existant
-          // //cout << "agregat deja existant" << endl;
-          if (!ancien) {
-            cto.insert(sw);
-          }
-          reinit_cycle(sw, tranobs);
-
-          sw.pop_back();
-          ancien = true;
-          // cout << "on enleve du chemin"<< endl;
-          delete reached_class;
-          e.first.first->Successors.insert(e.first.first->Successors.begin(),
-                                           Edge(pos, t));
-          pos->Predecessors.insert(pos->Predecessors.begin(),
-                                   Edge(e.first.first, t));
-          // //cout<<"les pred existe sont avec "<<Edge(e.first.first,
-          // t).second+1<<endl;
-          g.addArc();
-        }
-      }
-
-    } else {
-      sw.pop_back();
-      ancien = true;
-    }
-  }
-
-  return cto;
-}
-///////////////////////////////////////////////
-
-/////////////////////////////////////////////
-void RdPBDD::complete_sog(MDGraph &g, map<int, int> tranobs) {
-  set<chem> cto;
-  vector<int> sw;
-  bool ancien = true;
-  set<int> cov;
-  typedef pair<Class_Of_State *, bdd> couple;
-  typedef pair<couple, Set> Pair;
-  typedef stack<Pair> pile;
-  pile st;
-  int t;
-
-  double d, tps;
-  d = (double)clock() / (double)CLOCKS_PER_SEC;
-  TabMeta = new bdd[100000];
-  nbmetastate = 0;
-  MaxIntBdd = 0;
-  NbIt = 0;
-  itext = itint = 0;
-  Class_Of_State *reached_class;
-  Set fire;
-
-  Class_Of_State *c = new Class_Of_State;  // construction duu premier agregat
-  {
-    bdd Complete_meta_state = Accessible_epsilon(M0);
-    fire = firable_obs(Complete_meta_state);
-    c->class_state = Complete_meta_state;
-    TabMeta[nbmetastate] = c->class_state;
-    nbmetastate++;
-    old_size = bdd_nodecount(c->class_state);
-    st.push(Pair(couple(c, Complete_meta_state), fire));
-  }
-  g.setInitialState(c);  // init
-  g.insert(c);
-  // cout << "premier agregat construit " << endl;
-  g.nbMarking += bdd_pathcount(c->class_state);
-  // for (auto f : fire) {
-  // cout << "le nb de fires de"<<c->class_state.id()<<" est " << fire.size()
-  // << " : t" << f+1 << endl;
-  // }
-  while (!st.empty()) {
-    // cout<<"il reste des elmts dans st " <<endl;
-    Pair e = st.top();
-    st.pop();
-    // nbmetastate--;
-
-    if (!e.second.empty())  // firable non vide
-    {
-      // //cout<<"il reste des fire" <<endl;
-      // choisir transition des transitions firables
-      int t = choisir(cov, e.second);
-      if (t != -1) {
-        ancien = false;
-        tranobs[t]--;
-        sw.push_back(t);
-        if (tranobs[t] == 0) {
-          cov.insert(t);
-          if (cov.size() == Observable.size()) {
-            if (!ancien) {
-              cto.insert(sw);
-            }
-            // return cto;
-          }
-        }
-      }
-
-      else {
-        t = *e.second.begin();
-        sw.push_back(t);
-      }
-      // cout<<"+++++++++++++++on a trouv� la transition t"<<t+1 <<" partant de
-      // " <<e.first.first->class_state.id()<<endl;
-      e.second.erase(t);
-
-      // //cout << "on ajoute t"<<t+1 <<" au chemin"<< endl;
-
-      // //cout << "ici on couvre  la transition t" << t+1 <<"qui a un second de
-      // "<<tranobs[t]<< " et le nb de trans couv = " << cov.size()
-      // <<"/"<<Observable.size()<< endl;
-      st.push(e);
-
-      double nbnode;
-      reached_class = new Class_Of_State;
-      {
-        // cout << "on traite un bdd avac la trans t" <<t+1<< endl;
-        bdd Complete_meta_state =
-            Accessible_epsilon(get_successor(e.first.second, t));
-        // cout << "on crre un bdd" << endl;
-        reached_class->class_state = Complete_meta_state;
-        Class_Of_State *pos = g.find(reached_class);
-        nbnode = bdd_pathcount(reached_class->class_state);
-
-        if (!pos)  // nouvel agregat
-        {
-          // cout << "nouvel agregat " << endl;
-          fire = firable_obs(Complete_meta_state);
-          // for (auto f : fire) {
-          // cout << "le nb de fires de  "<<reached_class->class_state.id()<<"
-          // est " << fire.size() << " : t" << f+1 << endl;
-          // }
-          // st.push(Pair(couple(reached_class, Complete_meta_state), fire));
-          nbmetastate++;
-          e.first.first->Successors.insert(e.first.first->Successors.begin(),
-                                           Edge(reached_class, t));
-          reached_class->Predecessors.insert(
-              reached_class->Predecessors.begin(), Edge(e.first.first, t));
-          // //cout<<"les pred sont avec t"<<Edge(e.first.first,
-          // t).second+1<<endl;
-          g.addArc();
-          g.insert(reached_class);  // ajouter l'agregat
-          if (fire.empty()) {
-            // cout << "pas de successeurs " << endl;
-
-            if (!ancien) {
-              cto.insert(sw);
-            }
-            reinit_cycle(sw, tranobs);
-
-            sw.pop_back();
-            // cout << "on enleve du chemin"<< endl;
-
-          } else {
-            // cout << "il y a  de successeurs " << endl;
-            st.push(Pair(couple(reached_class, Complete_meta_state), fire));
-          }
-
-        } else {  // agregat existant
-          // cout << "agregat"<< pos->class_state.id()<<" deja existant" <<
-          // endl;
-          if (!ancien) {
-            cto.insert(sw);
-          }
-          reinit_cycle(sw, tranobs);
-
-          sw.pop_back();
-          ancien = true;
-          // cout << "on enleve du chemin"<< endl;
-          delete reached_class;
-          e.first.first->Successors.insert(e.first.first->Successors.begin(),
-                                           Edge(pos, t));
-          pos->Predecessors.insert(pos->Predecessors.begin(),
-                                   Edge(e.first.first, t));
-          // cout<<"les pred existe sont avec "<<Edge(e.first.first,
-          // t).second+1<<endl;
-          g.addArc();
-        }
-      }
-
-    } else {
-      sw.pop_back();
-      ancien = true;
-    }
-  }
-}
-
-/////////////////////////////////////////////
-
-stack<pair<Class_Of_State *, bdd>> RdPBDD::recherche_points_entree(chem ch,
-                                                                   MDGraph &g) {
-  stack<pair<Class_Of_State *, bdd>> pt_entr;
-  pair<Class_Of_State *, bdd> p;
-  Class_Of_State *agr = new Class_Of_State;
-
-  bdd entree = M0;
-  agr = g.initialstate;
-
-  p.first = agr;
-  p.second = entree;
-  pt_entr.push(p);
-  // cout<<"l' 1 agr "<<agr->class_state.id()<<endl;
-  // cout<<"le 1 bdd "<<entree<<endl;
-
-  // //cout<<"les points de sortie avec  le t"<< ch[0]<<" sont :
-  // "<<FrontiereNodes1(agr->class_state,ch[0])<<endl;
-
-  for (std::vector<int>::iterator k = ch.begin(); k != ch.end() - 1; k++) {
-    int t = *k;
-    entree = relation[t](p.first->class_state);
-    p.second = entree;
-
-    for (auto succ : agr->Successors) {
-      if (succ.second == t) {
-        agr = succ.first;
-        // cout<<"avec trans "<<t+1<<endl;
-        // cout<<"l'id de agr "<<agr->class_state.id()<<endl;
-        break;
-      }
-    }
-
-    p.first = agr;
-    pt_entr.push(p);
-  }
-  // cout<<"entrey points found"<<endl;
-  return pt_entr;
-}
-
-////////////////////////////////////////////
-
-chem RdPBDD::chem_abs(chem ch, MDGraph &g) {
-  vector<int> ch_abstrait, ch_agregat;
-  int trans;
-  pair<Class_Of_State *, bdd> agr_entree;
-  stack<pair<Class_Of_State *, bdd>> point_entree;
-  bdd cible, source;
-  Class_Of_State *agr = new Class_Of_State;
-  // agr=dernier agr
-  point_entree = recherche_points_entree(ch, g);
-  int i = 0;
-  while (!point_entree.empty()) {
-    trans = *(ch.end() - 1);
-    // cout<<"on traite "<< trans+1<<endl;
-    agr_entree = point_entree.top();
-    // cout<<" l'agr trait� "<<agr_entree.first->class_state.id()<<endl;
-    point_entree.pop();  // effacer l'element le dernier element
-    cible = agr_entree.second;
-    agr = agr_entree.first;
-    if (i == 0) {
-      source = FrontiereNodes1(agr->class_state, trans);
-    }
-    if (i == 1) {
-      source = (relation[trans])[source];
-    }
-    ch_agregat = Sub_path_agregate(&source, cible, agr);
-    // cout<<"sub_path found "<<endl;
-    ch_abstrait.insert(ch_abstrait.begin(), trans);
-    ch_abstrait.insert(ch_abstrait.begin(), ch_agregat.begin(),
-                       ch_agregat.end());
-    ch.pop_back();
-    i = 1;
-  }
-
-  return ch_abstrait;
-}
-
-///////----------------- trouver un chemain de la source � cible dans un agregat
-/// agr---
-vector<int> RdPBDD::Sub_path_agregate(bdd *source, bdd cible,
-                                      Class_Of_State *agr) {
-  vector<int> ch_agregat;
-
-  pair<int, bdd> couple;
-  int transition;
-  int i = 0;
-  bdd courant = *source;
-  bdd egalite = cible & courant;
-
-  // cout<<"a difference between source and target "<< (egalite ==
-  // bdd_false())<<endl;
-  while (egalite == bdd_false()) {
-    couple = StepBackward1(courant, agr);
-    ch_agregat.insert(ch_agregat.begin(), couple.first);
-    courant = couple.second;
-    egalite = (cible & courant);
-
-    // cout<<"a difference between source and target  "<< (egalite ==
-    // bdd_false())<<endl;
-  }
-
-  *source = courant;
-  return ch_agregat;
-}
-
-/*-----------------------CanonizeR()----------------*/
-bdd RdPBDD::CanonizeR(bdd s, unsigned int i) {
-  bdd s1, s2;
-  do {
-    itext++;
-    s1 = s - bdd_nithvar(2 * i);
-    s2 = s - s1;
-    if ((!(s1 == bddfalse)) && (!(s2 == bddfalse))) {
-      bdd front = s1;
-      bdd reached = s1;
-      do {
-        // cout<<"premiere boucle interne \n";
-        itint++;
-        front = StepForward(front) - reached;
-        reached = reached | front;
-        s2 = s2 - front;
-      } while ((!(front == bddfalse)) && (!(s2 == bddfalse)));
-    }
-    if ((!(s1 == bddfalse)) && (!(s2 == bddfalse))) {
-      bdd front = s2;
-      bdd reached = s2;
-      do {
-        // //cout<<"deuxieme boucle interne \n";
-        itint++;
-        front = StepForward(front) - reached;
-        reached = reached | front;
-        s1 = s1 - front;
-      } while ((!(front == bddfalse)) && (!(s1 == bddfalse)));
-    }
-    s = s1 | s2;
-    i++;
-  } while ((i < Nb_places) && ((s1 == bddfalse) || (s2 == bddfalse)));
-  if (i >= Nb_places) {
-    // cout<<"____________oooooooppppppppsssssssss________\n";
-    return s;
-  } else {
-    // cout<<"________________p a s o o o p p p s s s ______\n";
-    return (CanonizeR(s1, i) | CanonizeR(s2, i));
-  }
-}
-/*---------------------------  Set_Bloc()  -------*/
-bool RdPBDD::Set_Bloc(bdd &S) const {
-  // cout<<"Ici detect blocage \n";
-  int k = 0;
-  bdd Mt = bddtrue;
-  for (vector<Trans>::const_iterator i = relation.begin(); i != relation.end();
-       i++, k++) {
-    Mt = Mt & !((*i).Precond);
-  }
-  return ((S & Mt) != bddfalse);
-  // BLOCAGE
-}
-/*-------------------------Set_Div() � revoir -----------------------------*/
-bool RdPBDD::Set_Div(bdd &S) const {
-  Set::const_iterator i;
-  bdd To, Reached;
-  // cout<<"Ici detect divergence \n";
-  Reached = S;
-  do {
-    bdd From = Reached;
-    for (i = NonObservable.begin(); !(i == NonObservable.end()); i++) {
-      To = relation[*i](Reached);
-      Reached = Reached | To;  // Reached=To ???
-      // Reached=To;
-    }
-    if (Reached == From)
-      // cout<<"SEQUENCE DIVERGENTE \n";
-      return true;
-    // From=Reached;
-  } while (Reached != bddfalse);
-  return false;
-  // cout<<"PAS DE SEQUENCE DIVERGENTE \n";
-}
-/*-----------FrontiereNodes() pour bdd ---------*/
-bdd RdPBDD::FrontiereNodes(bdd From) const {
-  bdd res = bddfalse;
-  for (Set::const_iterator i = Observable.begin(); !(i == Observable.end());
-       i++)
-    res = res | (From & relation[*i].Precond);
-  return res;
-}
-/*-----------FrontiereNodes1() pour bdd ---------*/
-bdd RdPBDD::FrontiereNodes1(bdd From, int t) {
-  bdd res = bddfalse;
-  res = res | (From & relation[t].Precond);
-  return res;
-}
-/*-------- Produit synchronis� � la vol�e de n graphes d'observation :
- * Adaptation � la capture des s�quences bloquantes et les s�quences
- * divergentes----------------------*/
-void RdPBDD::GeneralizedSynchProduct1(Modular_Obs_Graph &Gv, int NbSubnets,
-                                      RdPBDD *Subnets[], int nbbddvar) {
-  // cout<<"_____________  GeneralizedSynchProduct1
-  // _________________________"<<NbSubnets<<"sous-reseaux "<<endl; for(int
-  // k=0;k<NbSubnets;k++)
-  //	//cout<<*Subnets[k]<<endl;
-  int pos_trans(TRANSITIONS, string);
-  TabMeta = new bdd[1000000];
-  nbmetastate = 0;
-  MaxIntBdd = 0;
-  nbmetastate = 0;
-  Stack st;
-  // int cpt=0;
-  int k;
-  bdd *Complete_meta_state = new bdd[NbSubnets];
-  Set *fire = new Set[NbSubnets];
-  Modular_Class_Of_State *Meta_State = new Modular_Class_Of_State;
-  for (k = 0; k < NbSubnets; k++) {
-    Complete_meta_state[k] = Subnets[k]->Accessible_epsilon(Subnets[k]->M0);
-    fire[k] = Subnets[k]->firable_obs(Complete_meta_state[k]);
-    // Meta_State->State.push_back(Subnets[k]->FrontiereNodes(Complete_meta_state[k]));
-    // Meta_State->State.push_back(Subnets[k]->CanonizeR(Subnets[k]->FrontiereNodes(Complete_meta_state[k]),0));
-    Meta_State->State.push_back(
-        Subnets[k]->CanonizeR(Complete_meta_state[k], 0));
-    /*-------------------- STAT ----------------------*/
-    TabMeta[nbmetastate] = Meta_State->State[k];
-    nbmetastate++;
-  }
-  old_size = bdd_anodecount(TabMeta, nbmetastate);
-  Meta_State->blocage = true;
-  for (k = 0; ((k < NbSubnets) && (Meta_State->blocage)); k++)
-    Meta_State->blocage =
-        Meta_State->blocage && Subnets[k]->Set_Bloc(Complete_meta_state[k]);
-  Meta_State->boucle = false;
-  for (k = 0; ((k < NbSubnets) && (!Meta_State->boucle)); k++)
-    Meta_State->boucle =
-        Meta_State->boucle || Subnets[k]->Set_Div(Complete_meta_state[k]);
-  Gv.setInitialState(Meta_State);
-  Gv.insert(Meta_State);
-  nbmetastate++;
-  st.push(StackElt(Couple(Meta_State, Complete_meta_state), fire));
-  do {
-    NbIt++;
-    StackElt e = st.top();
-    st.pop();
-    for (k = 0; k < NbSubnets; k++) {
-      while (!e.second[k].empty()) {
-        int t = *e.second[k].begin();
-        e.second[k].erase(t);
-        bool ok = true;
-        Set ConcernedSubnets;
-        ConcernedSubnets.insert(k);
-        string tmp = Subnets[k]->transitions[t].name;
-        for (int j = 0; j < NbSubnets; j++) {
-          if (j != k) {
-            int num = Subnets[j]->transitionName[tmp];
-            int pos = pos_trans(Subnets[j]->transitions, tmp);
-            if ((pos != -1) && !(Subnets[j]->Observable.find(num) ==
-                                 Subnets[j]->Observable.end())) {
-              ConcernedSubnets.insert(j);
-              if (e.second[j].find(num) == e.second[j].end()) {
-                ok = false;
-              } else
-                e.second[j].erase(num);
-            }
-          }
-        }
-        if (ok) {
-          Complete_meta_state = new bdd[NbSubnets];
-          fire = new Set[NbSubnets];
-          Meta_State = new Modular_Class_Of_State;
-          for (int j = 0; j < NbSubnets; j++) {
-            if (ConcernedSubnets.find(j) == ConcernedSubnets.end()) {
-              Complete_meta_state[j] = e.first.second[j];
-              Meta_State->State.push_back(e.first.first->State[j]);
-            } else {
-              Complete_meta_state[j] =
-                  Subnets[j]->Accessible_epsilon(Subnets[j]->get_successor(
-                      e.first.second[j], Subnets[j]->transitionName[tmp]));
-              // Point de sortie
-              // Meta_State->State.push_back(Subnets[j]->FrontiereNodes(Complete_meta_state[j]));
-              // Meta_State->State.push_back(Subnets[j]->CanonizeR(Subnets[j]->FrontiereNodes(Complete_meta_state[j]),0));
-              Meta_State->State.push_back(
-                  Subnets[j]->CanonizeR(Complete_meta_state[j], 0));
-              /*-------------------- STAT ----------------------*/
-              TabMeta[nbmetastate] = Meta_State->State[k];
-              nbmetastate++;
-            }
-            fire[j] = Subnets[j]->firable_obs(Complete_meta_state[j]);
-          }
-          Modular_Class_Of_State *pos = Gv.find(Meta_State);
-          if (!pos) {
-            old_size = bdd_anodecount(TabMeta, nbmetastate);
-            // Calcul de deadlock et loop attributes
-            Meta_State->blocage = true;
-            for (int j = 0; ((j < NbSubnets) && (Meta_State->blocage)); j++)
-              Meta_State->blocage &=
-                  Subnets[j]->Set_Bloc(Complete_meta_state[j]);
-            Meta_State->boucle = false;
-            for (int j = 0; ((j < NbSubnets) && (!Meta_State->boucle)); j++)
-              Meta_State->boucle |= Subnets[j]->Set_Div(Complete_meta_state[j]);
-            st.push(StackElt(Couple(Meta_State, Complete_meta_state), fire));
-            e.first.first->Successors.insert(e.first.first->Successors.begin(),
-                                             Modular_Edge(Meta_State, tmp));
-            Meta_State->Predecessors.insert(Meta_State->Predecessors.begin(),
-                                            Modular_Edge(e.first.first, tmp));
-            Gv.addArc();
-            Gv.insert(Meta_State);
-          } else {
-            e.first.first->Successors.insert(e.first.first->Successors.begin(),
-                                             Modular_Edge(pos, tmp));
-            pos->Predecessors.insert(pos->Predecessors.begin(),
-                                     Modular_Edge(e.first.first, tmp));
-            Gv.addArc();
-            delete Meta_State;
-            // Neoud d�ja rencontr� ;
-          }
-        }
-      }
-    }
-  } while (!st.empty());
-  // cout<<" MAXIMAL INTERMEDIARY BDD SIZE \n"<<MaxIntBdd<<endl;
-  // cout<<"OLD SIZE : "<<bdd_anodecount(TabMeta,nbmetastate)<<endl;
-  // cout<<"NB SHARED NODES : "<<bdd_anodecount(TabMeta,nbmetastate)<<endl;
-  // cout<<"NB META STATE DANS CONSTRUCTION : "<<nbmetastate<<endl;
-  // cout<<"NB ITERATIONS CONSTRUCTION : "<<NbIt<<endl;
-  // cout<<"NB ITERATIONS EXTERNES : "<<itext<<endl;
-  // cout<<"NB ITERATIONS INTERNES : "<<itint<<endl;
-}
-/*------------------------EmersonLey ----------------------------*/
-bdd RdPBDD::EmersonLey(bdd S, bool trace) {
-  // cout<<"ICI EMERSON LEY \n";
-  double TpsInit, TpsDetect;
-  double debitext, finitext;
-  TpsInit = (double)(clock()) / CLOCKS_PER_SEC;
-  bdd b = S;
-  bdd Fair = bdd_ithvar(2 * Nb_places - 1);
-  // cout<<"PLACE TEMOIN \n";
-  // cout<<places[places.size()-1].name<<endl;
-  bdd oldb;
-  bdd oldd, d;
-  int extit = 0;
-  int init = 0;
-  do {
-    extit++;
-    if (trace) {
-      // cout<<"ITERATION EXTERNES NUMERO :"<<extit<<endl;
-      debitext = (double)(clock()) / CLOCKS_PER_SEC;
-      // cout<<"TAILLE DE B AVANT IT INTERNE : "<<bdd_nodecount(b)<<endl;
-      // cout<<endl<<endl;
-    }
-    oldb = b;
-    // cout<<"Fair : "<<Fair.id()<<endl;
-    d = b & Fair;
-    // cout<<"d : "<<d.id()<<endl;
-    // init=0;
-    do {
-      init++;
-      if (trace) {
-        // cout<<"ITERATION INTERNES NUMERO :"<<init<<endl;
-        // cout<<"HEURE : "<<(double)(clock()) / CLOCKS_PER_SEC<<endl;
-        // cout<<"TAILLE DE D : "<<bdd_nodecount(d)<<endl;
-      }
-      oldd = d;
-      bdd inter = b & StepForward2(d);
-      // cout<<"Tille de inter :"<<bdd_nodecount(inter)<<endl;
-      d = d | inter;
-    } while (!(oldd == d));
-    if (trace)
-      // cout<<"\nTAILLE DE D APRES ITs INTERNES : "<<bdd_nodecount(d)<<endl;
-      b = b & StepBackward2(d);
-    init++;
-    if (trace) {
-      // cout<<"\n\nTAILLE DE B APRES ELEMINER LES PRED DE D :
-      // "<<bdd_nodecount(b)<<endl;
-      finitext = (double)(clock()) / CLOCKS_PER_SEC;
-      // cout<<"DUREE DE L'ITERATION EXTERNE NUMERO "<<extit<<"  :
-      // "<<finitext-debitext<<endl;
-      // cout<<endl<<"_________________________________________________\n\n";
-    }
-  } while (!(b == oldb));
-  // cout<<"NOMBRE D'ITERATIONS EXTERNES -----:"<<extit<<endl;
-  // cout<<"NOMBRE D'ITERATIONS INTERNES -----:"<<init<<endl;
-  TpsDetect = ((double)(clock()) / CLOCKS_PER_SEC) - TpsInit;
-  // cout << "DETECTION DE CYCLE TIME  " << TpsDetect << endl;
-  return b;
-}
diff --git a/src/RdPBDD.hpp b/src/RdPBDD.hpp
deleted file mode 100644
index 999a9874782739a3826830825cd5d7c91dfc6639..0000000000000000000000000000000000000000
--- a/src/RdPBDD.hpp
+++ /dev/null
@@ -1,95 +0,0 @@
-/* -*- C++ -*- */
-#ifndef RdPBDD_H
-#define RdPBDD_H
-#include <stack>
-#include <vector>
-
-#include "MDGraph.hpp"
-#include "Modular_Class_of_state.hpp"
-#include "Modular_Obs_Graph.hpp"
-#include "Net.hpp"
-#include "bdd.h"
-// #include"Modular_Obs_Graph.hpp"
-/***********************************************/
-
-class Trans {
- public:
-  Trans(const bdd& var, bddPair* pair, const bdd& rel, const bdd& prerel,
-        const bdd& Precond, const bdd& Postcond);
-  bdd operator()(const bdd& op) const;
-  bdd operator[](const bdd& op) const;
-  friend class RdPBDD;
-
- private:
-  bdd var;
-  bddPair* pair;
-  bdd Precond;
-  bdd Postcond;
-  bdd rel;
-  bdd prerel;
-};
-typedef set<int> Set;
-typedef vector<int> chem;
-
-class RdPBDD {
- private:
-  vector<class Transition> transitions;
-  Set Observable;
-  Set NonObservable;
-  map<string, int> transitionName;
-  unsigned int Nb_places;
-
- public:
-  bdd M0;
-  bdd currentvar;
-  vector<Trans> relation;
-  bdd ReachableBDD1();
-  bdd ReachableBDD2();
-  bdd ReachableBDD3();
-  /* G�n�ration de graphes d'observation */
-  stack<pair<Class_Of_State*, bdd>> recherche_points_entree(chem ch,
-                                                            MDGraph& g);
-  set<chem> chem_obs(MDGraph& g, map<int, int> tranobs);
-  void complete_sog(MDGraph& g, map<int, int> tranobs);
-  chem chem_abs(chem ch, MDGraph& g);
-  vector<int> Sub_path_agregate(bdd* source, bdd cible, Class_Of_State* agr);
-  void bdd_firable_obs(Class_Of_State* agr, int t);
-
-  void compute_canonical_deterministic_graph_Opt(MDGraph& g);
-  bdd Accessible_epsilon(bdd From);
-  bdd Accessible_epsilon2(bdd From);
-  bdd StepForward(bdd From);
-  bdd StepForward2(bdd From);
-  bdd StepBackward(bdd From);
-  pair<int, bdd> StepBackward1(bdd From, Class_Of_State* agr);  // ajout�
-  bdd StepBackward2(bdd From);
-  bdd EmersonLey(bdd S, bool trace);
-  Set firable_obs(bdd State);
-  bdd get_successor(bdd From, int t);
-  void GeneralizedSynchProduct1(Modular_Obs_Graph& Gv, int NbSubnets,
-                                RdPBDD* Subnets[], int nbbddvar);
-  bool Set_Div(bdd& S) const;
-  bool Set_Bloc(bdd& S) const;
-  bdd FrontiereNodes(bdd From) const;
-  bdd FrontiereNodes1(bdd From, int t);  // ajout�
-  bdd CanonizeR(bdd s, unsigned int i);
-  RdPBDD(const net&, map<int, int> observables, Set NonObservables,
-         int BOUND = 32, bool init = false);
-
-  ~RdPBDD(){};
-};
-/*____________Structure utiles aux produit synchronis� g�n�ralis� de n graphes
- * d'observations ________*/
-
-typedef pair<Modular_Class_Of_State*, bdd*> Couple;
-typedef pair<Couple, Set*> StackElt;
-typedef stack<StackElt> Stack;
-typedef vector<StackElt> PileInt;
-typedef pair<StackElt, PileInt> CanonicRep;
-// typedef stack<CanonicRep> STACK;
-typedef stack<PileInt> STACK;
-int InStack(PileInt, Modular_Class_Of_State*);
-bool MemeVecteur(vector<bdd>, vector<bdd>);
-bool Inferieur(vector<bdd>, vector<bdd>);
-
-#endif
diff --git a/src/aggregate.hpp b/src/aggregate.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..f2695b9345f6dcca8b246097b60e6d125071f205
--- /dev/null
+++ b/src/aggregate.hpp
@@ -0,0 +1,27 @@
+#ifndef AGGREGATE_H_
+#define AGGREGATE_H_
+
+#include <set>
+#include <vector>
+
+#include "bdd.h"
+
+// definition of types
+class Aggregate;
+typedef std::pair<Aggregate*, int> Edge;
+typedef std::vector<Edge> Edges;
+
+class Aggregate {
+ public:
+  Aggregate() = default;
+
+  bdd state;
+  bool boucle{false};
+  bool blocage{false};
+  bool visited{false};
+
+  std::vector<Edge> predecessors;
+  std::vector<Edge> successors;
+};
+
+#endif  // AGGREGATE_H_
diff --git a/src/main.cpp b/src/main.cpp
index 73bbdd2c63b1464bf7e532c42394400b89fa37a9..7d4bc844fca9f1419210a944065385a110688ca0 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -1,8 +1,5 @@
-#include <time.h>
-
 #include <CLI11.hpp>
-#include <algorithm>
-#include <cmath>
+#include <ctime>
 #include <fstream>
 #include <iomanip>
 #include <iostream>
@@ -10,147 +7,149 @@
 #include <string>
 
 #include "Net.hpp"
-#include "RdPBDD.hpp"
 #include "bdd.h"
-#include "fdd.h"
+#include "petri_net_sog.hpp"
 
-using namespace std;
+// BDD initial values
+constexpr int BDD_INITIAL_NUM_NODES_ = 10000;
+constexpr int BDD_SIZE_CACHES = 10000;
 
 // definition of type path
-typedef set<vector<int>> paths_t;
+typedef std::vector<int> path_t;
+typedef std::set<path_t> paths_t;
 
 /**
- * @brief Return the current time in secods
+ * Return the current time in seconds
  * @return time in seconds
  */
-double getTime() { return (double)clock() / (double)CLOCKS_PER_SEC; }
+double GetTime() {
+  return static_cast<double>(clock()) / static_cast<double>(CLOCKS_PER_SEC);
+}
 
 /**
- * @brief Save in a file the information about observable paths
+ * Save in a file the information about observable paths
  * @param output_file Path to the file where the output will be saved
  * @param obs_paths Set of observable paths
  */
-void save_observable_paths(const paths_t& obs_paths,
-                           const string& output_file) {
-  cout << "\nSaving results in  " << output_file << endl;
+void SaveObservablePaths(const paths_t& obs_paths, const string& output_file) {
+  cout << "\nSaving results in  " << output_file << '\n';
 
   ofstream my_file(output_file);
-  my_file << "# paths: " << obs_paths.size() << endl;
+  my_file << "# paths: " << obs_paths.size() << '\n';
 
-  int n = 1;
-  for (auto path : obs_paths) {
-    my_file << "path #" << n << ": " << path.size() << " transitions" << endl;
-    n++;
+  int path_id = 1;
+  for (const vector<int>& path : obs_paths) {
+    my_file << "path #" << path_id << ": " << path.size() << " transitions"
+            << '\n';
+    path_id++;
   }
 }
 
 /**
- * @brief Print statistics about the generated abstract paths
+ * Print statistics about the generated abstract paths
  * @param abstract_paths Set of abstract paths
  */
-void print_stats(const paths_t& abstract_paths) {
-  int sum_transtions = 0;
+void PrintStats(const paths_t& abstract_paths) {
+  size_t sum_transitions = 0;
   set<int> transitions;
 
-  for (auto path : abstract_paths) {
-    for (auto tr : path) {
-      transitions.insert(tr);
+  for (const path_t& path : abstract_paths) {
+    for (int transition_id : path) {
+      transitions.insert(transition_id);
     }
-    sum_transtions += path.size();
+    sum_transitions += path.size();
   }
 
-  float nb_paths = abstract_paths.size();
-  float average = sum_transtions / nb_paths;
+  const auto nb_paths = static_cast<double>(abstract_paths.size());
+  const double average = static_cast<double>(sum_transitions) / nb_paths;
 
-  float std_deviation = 0;
-  for (auto path : abstract_paths) {
+  double std_deviation = 0;
+  for (const path_t& path : abstract_paths) {
     std_deviation += pow(path.size() - average, 2);
   }
   std_deviation = sqrt(std_deviation / nb_paths);
 
-  cout << "# abstract paths: " << nb_paths << endl;
-  cout << "average # of transitions per abstract path: " << average << endl;
-  cout << "standard deviation: " << std_deviation << endl;
-  cout << "# covered transitions: " << transitions.size() << endl;
+  cout << "# abstract paths: " << nb_paths << '\n';
+  cout << "average # of transitions per abstract path: " << average << '\n';
+  cout << "standard deviation: " << std_deviation << '\n';
+  cout << "# covered transitions: " << transitions.size() << '\n';
 }
 
 /**
- * @brief Prints the output of the tool
+ * Prints the output of the tool
  * @param model Petri net model
- * @param g SOG graph
+ * @param sog SOG graph
  * @param obs_transitions observable transitions
  * @param obs_paths set of observable paths
  * @param abstract_paths set of abstract paths
- * @param elapsed_time time consumed by the computation
  * @param output_file file where the information will be saved
  */
-void print_output(const net& model, MDGraph& g,
-                  const map<int, int>& obs_transitions,
-                  const paths_t& obs_paths, const paths_t& abstract_paths,
-                  const string& output_file) {
+void PrintOutput(const net& model, SOG& sog,
+                 const map<int, int>& obs_transitions, const paths_t& obs_paths,
+                 const paths_t& abstract_paths, const string& output_file) {
   // print SOG information
-  g.printCompleteInformation();
-  cout << "\n# transition: " << model.transitions.size() << endl;
-  cout << "# places: " << model.places.size() << endl;
-  cout << "# observable transitions: " << obs_transitions.size() << endl
-       << endl;
+  sog.PrintCompleteInformation();
+  cout << "\n# transition: " << model.transitions.size() << '\n';
+  cout << "# places: " << model.places.size() << '\n';
+  cout << "# observable transitions: " << obs_transitions.size() << "\n\n";
 
   // print stats
-  print_stats(abstract_paths);
+  PrintStats(abstract_paths);
 
   // save observable paths in the folder
-  save_observable_paths(obs_paths, output_file);
+  SaveObservablePaths(obs_paths, output_file);
 }
 
 /**
- * @brief Returns the name of the model from the filename
+ * Returns the name of the model from the filename
  * @param filename Path to the petri net model
  * @return name of the model
  */
-string get_model_name(const string& filename) {
-  int pp = filename.find(".n") + 1;
-  int ps = filename.rfind("/") + 1;
-  string name = filename.substr(ps, pp - ps - 1);
+string GetModelName(const string& filename) {
+  const auto end_pos = filename.find(".n") + 1;
+  const auto start_pos = filename.rfind('/') + 1;
+  string name = filename.substr(start_pos, end_pos - start_pos - 1);
   return name;
 }
 
 /**
- * @brief Load a petri net from a .net file
+ * Load a petri net from a .net file
  * @param filename Path to the petri net model
  * @return Petri net model
  */
-net load_net(const string& filename) {
+net LoadPetriNet(const string& filename) {
   cout << "Parsing net: " << filename << " ... ";
   net model(filename.c_str());
-  cout << "done" << endl;
+  cout << "done\n";
 
   return model;
 }
 
 /**
- * @brief Load observable transitions from a file
+ * Load observable transitions from a file
  * @param model Petri net model
  * @param file Path to the file with the observable transitions
  * @param obs_trans set of observable transitions
  * @param unobs_trans set of unobservable transitions
  */
-void load_transitions(const net& model, const string& file,
-                      map<int, int>& obs_trans, set<int>& unobs_trans) {
-  string line;
+void LoadTransitions(const net& model, const string& file,
+                     map<int, int>& obs_trans, set<int>& unobs_trans) {
   ifstream my_file(file);
 
   // TODO: What means the second element of obs_trans type ? {<tr_id>, 1}
   if (my_file) {
+    string line;
     while (getline(my_file, line)) {
       // TODO: catch error when the file has transition different to t
-      int tr = stoi(line.substr(1));  // for a transition t1, we take the id 1
-      obs_trans.insert({tr, 1});
+      // for a transition t1, we take the id 1
+      int trans_id = stoi(line.substr(1));
+      obs_trans.insert({trans_id, 1});
     }
   }
 
   // TODO: What happens if the transition started in t1 instead of t0 ?
   // compute the unobservable transitions
-  for (long unsigned int i = 0; i < model.transitions.size(); i++) {
+  for (auto i = 0; i < model.transitions.size(); i++) {
     if ((obs_trans.find(i)) == obs_trans.end()) {
       unobs_trans.insert(i);
     }
@@ -158,19 +157,19 @@ void load_transitions(const net& model, const string& file,
 }
 
 /**
- * @brief Find the observable transitions needed for covering all the model
+ * Find the observable transitions needed for covering all the model
  * @param model Petri net model
  * @param obs_trans set of observable transitions
  * @param unobs_trans set of unobservable transitions
  */
-void find_observable_transitions(net& model, map<int, int>& obs_trans,
-                                 set<int>& unobs_trans) {
+void FindObservableTransitions(net& model, map<int, int>& obs_trans,
+                               set<int>& unobs_trans) {
   // compute the unobservable transitions of the model using the pattern
   // TODO: why cannot be const model ?
   unobs_trans = model.calcul1();
 
   // compute the observable transitions
-  for (long unsigned int i = 0; i < model.transitions.size(); i++) {
+  for (auto i = 0; i < model.transitions.size(); i++) {
     if ((unobs_trans.find(i)) == unobs_trans.end()) {
       obs_trans.insert({i, 1});
     }
@@ -178,90 +177,93 @@ void find_observable_transitions(net& model, map<int, int>& obs_trans,
 }
 
 /**
- * @brief Generate abstract paths of a Petri net
+ * Generate abstract paths of a Petri net
  * @param net_file Path to file of the petri net model
  * @param bound SOG bound
  * @param transitions_file Path to file of observable transitions
  * @param output_folder Path to folder where output files will be saved
  */
-void compute_abstract_paths(const string& net_file, int bound,
-                            const string& transitions_file,
-                            const string& output_folder) {
-  MDGraph g;
+void ComputeAbstractPaths(const string& net_file, int bound,
+                          const string& transitions_file,
+                          const string& output_folder) {
+  SOG sog;
   paths_t obs_paths;
   paths_t abstract_paths;
   set<int> unobs_trans;
   map<int, int> obs_trans;
-  net model = load_net(net_file);
+  net model = LoadPetriNet(net_file);
+
+  // BDD initialization.
+  // See https://buddy.sourceforge.net/manual/group__kernel.html
+  bdd_init(BDD_INITIAL_NUM_NODES_, BDD_SIZE_CACHES);
 
-  // BDD initialization
-  bdd_init(10000, 10000);
+  // Suppress GC messages
+  bdd_gbc_hook(nullptr);
 
-  float start_time = getTime();
+  auto start_time = GetTime();
 
   // if a path with transitions is not given, then we apply the algorithm to
   // find the needed observable transitions to cover all the behaviors
-  if (transitions_file == "") {
+  if (transitions_file.empty()) {
     cout << "\nComputing observable transitions ...";
-    float start_obs_time = getTime();
-    find_observable_transitions(model, obs_trans, unobs_trans);
-    float obs_time = getTime() - start_obs_time;
+    auto start_obs_time = GetTime();
+    FindObservableTransitions(model, obs_trans, unobs_trans);
+    auto obs_time = GetTime() - start_obs_time;
     cout << " done\nTime for computing observable transitions: " << obs_time
-         << " seconds" << endl;
+         << " seconds\n";
   } else {
-    load_transitions(model, transitions_file, obs_trans, unobs_trans);
+    LoadTransitions(model, transitions_file, obs_trans, unobs_trans);
   }
 
   // build the net
   cout << "\nBuilding net ...";
-  float start_net_time = getTime();
-  RdPBDD R(model, obs_trans, unobs_trans, bound, true);
-  float net_time = getTime() - start_net_time;
-  cout << " done\nTime for computing the net: " << net_time << " seconds"
-       << endl;
+  auto start_net_time = GetTime();
+  PetriNetSOG pn_sog(model, obs_trans, unobs_trans, bound, true);
+  auto net_time = GetTime() - start_net_time;
+  cout << " done\nTime for computing the net: " << net_time << " seconds\n";
 
   // compute the observable paths
   cout << "\nComputing observable paths ...";
-  float start_paths_time = getTime();
-  obs_paths = R.chem_obs(g, obs_trans);
-  float paths_time = getTime() - start_paths_time;
+  auto start_paths_time = GetTime();
+  obs_paths = pn_sog.ObservablePaths(sog, obs_trans);
+  auto paths_time = GetTime() - start_paths_time;
   cout << " done\nTime for computing observable paths: " << paths_time
-       << " seconds" << endl;
+       << " seconds\n";
 
   // add abstract paths
   cout << "\nComputing abstract paths ...";
-  float start_abstract_time = getTime();
-  for (auto path : obs_paths) {
-    abstract_paths.insert(R.chem_abs(path, g));
+  auto start_abstract_time = GetTime();
+  for (const path_t& path : obs_paths) {
+    abstract_paths.insert(pn_sog.AbstractPaths(path, sog));
   }
-  float abstract_time = getTime() - start_abstract_time;
+  auto abstract_time = GetTime() - start_abstract_time;
   cout << " done\nTime for computing abstract paths: " << abstract_time
-       << " seconds" << endl;
+       << " seconds\n";
 
   // time for generating the paths
-  float elapsed_time = getTime() - start_time;
-  cout << "\nTotal time: " << elapsed_time << " seconds" << endl;
+  auto elapsed_time = GetTime() - start_time;
+  cout << "\nTotal time: " << elapsed_time << " seconds\n";
 
   // print output
-  string model_name = get_model_name(net_file);
+  string model_name = GetModelName(net_file);
   string output_file = output_folder + "/obs_paths_" + model_name + ".txt";
-  print_output(model, g, obs_trans, obs_paths, abstract_paths, output_file);
+  PrintOutput(model, sog, obs_trans, obs_paths, abstract_paths, output_file);
 }
 
 /******************************************************************************
  * Main function
  ******************************************************************************/
-int main(int argc, char** argv) {
+int main(const int argc, char** argv) {
   CLI::App app{
       "sogMBT: Symbolic Observation Graph-Based Generation of Test Paths"};
 
-  string input_file = "";
+  string input_file;
   app.add_option("--input-net", input_file, "Petri net file")
       ->type_name("Path")
       ->required()
       ->check(CLI::ExistingFile);
 
-  string output_folder = "";
+  string output_folder;
   app.add_option("--output-folder", output_folder, "output folder")
       ->type_name("Path")
       ->required()
@@ -281,7 +283,7 @@ int main(int argc, char** argv) {
   // parse arguments
   CLI11_PARSE(app, argc, argv);
 
-  compute_abstract_paths(input_file, bound, obs_file, output_folder);
+  ComputeAbstractPaths(input_file, bound, obs_file, output_folder);
 
   return 0;
 }
diff --git a/src/modular_aggregate.hpp b/src/modular_aggregate.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..bb0100fb439a9be0e461d859e6f6bf48bae36d7c
--- /dev/null
+++ b/src/modular_aggregate.hpp
@@ -0,0 +1,49 @@
+#ifndef MODULAR_AGGREGATE_H_
+#define MODULAR_AGGREGATE_H_
+
+#include <set>
+#include <vector>
+
+#include "bdd.h"
+
+class ModularAggregate;
+typedef std::pair<ModularAggregate *, std::string> ModularEdge;
+
+// TODO: check this
+struct LessModularEdge {
+  bool operator()(const ModularEdge &a1, const ModularEdge &a2) const {
+    return (a1.first < a2.first);
+  }
+};
+
+typedef std::set<ModularEdge, LessModularEdge> ModularEdges;
+
+class ModularAggregate {
+ public:
+  ModularAggregate() = default;
+
+  std::vector<bdd> state;
+  bool boucle{false};
+  bool blocage{false};
+  bool visited{false};
+
+  std::set<ModularEdge> predecessors;
+  std::set<ModularEdge> successors;
+
+  /**
+   * Operator << overload
+   * @param os
+   * @param c
+   * @return
+   */
+  friend std::ostream &operator<<(std::ostream &os, const ModularAggregate &c) {
+    os << "{";
+    for (const auto &k : c.state) {
+      os << k.id() << ", ";
+    }
+    os << "}";
+    return (os);
+  }
+};
+
+#endif  // MODULAR_AGGREGATE_H_
diff --git a/src/modular_sog.cpp b/src/modular_sog.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..1f2d3fa01458d1358ca841f22a3ce66ddb41ed69
--- /dev/null
+++ b/src/modular_sog.cpp
@@ -0,0 +1,173 @@
+#include "modular_sog.hpp"
+
+// intermediate table to calculate the size (nb bdd) of the graph
+bdd *tmp_table;
+
+void ModularSOG::set_initial_state(ModularAggregate *s) {
+  current_state = initial_state = s;
+}
+
+ModularAggregate *ModularSOG::FindState(const ModularAggregate *s) const {
+  for (auto *GONode : GONodes) {
+    bool arret = false;
+    for (auto k = 0; ((k < (s->state).size()) && (!arret)); k++) {
+      if ((s->state[k] == GONode->state[k]) == 0) {
+        arret = true;
+      }
+    }
+
+    if (!arret) {
+      return GONode;
+    }
+  }
+
+  return nullptr;
+}
+
+ModularAggregate *ModularSOG::FindState2(const ModularAggregate *s,
+                                            Set states_involved) const {
+  for (auto *GONode : GONodes) {
+    bool arret = false;
+    for (auto k = 0; ((k < (s->state).size()) && (!arret)); k++) {
+      if (states_involved.find(k) != states_involved.end()) {
+        if ((s->state[k] == GONode->state[k]) == 0) {
+          arret = true;
+        }
+      }
+    }
+
+    if (!arret && (s->blocage == GONode->blocage) &&
+        (s->boucle == GONode->boucle)) {
+      return GONode;
+    }
+  }
+
+  return nullptr;
+}
+
+void ModularSOG::AddArc(ModularAggregate *pred, ModularAggregate *succ,
+                        const char *trans) {
+  const auto arc = ModularEdge(succ, trans);
+  const auto cra = ModularEdge(pred, trans);
+  if (pred->successors.find(arc) == pred->successors.end()) {
+    pred->successors.insert(pred->successors.begin(), arc);
+    succ->predecessors.insert(succ->predecessors.begin(), cra);
+    nb_arcs++;
+  }
+}
+
+void ModularSOG::InsertState(ModularAggregate *s) {
+  s->visited = false;
+  this->GONodes.push_back(s);
+  nb_states++;
+}
+
+void ModularSOG::InitVisit(ModularAggregate *s, size_t nb) {
+  if (nb <= nb_states) {
+    s->visited = false;
+    for (const auto &successor : s->successors) {
+      if (successor.first->visited) {
+        nb++;
+        InitVisit(successor.first, nb);
+      }
+    }
+  }
+}
+
+void ModularSOG::TabBDDNodes(ModularAggregate *s, size_t &nb) {
+  if (!s->visited) {
+    for (unsigned int k = 0; k < s->state.size(); k++, nb++) {
+      tmp_table[nb - 1] = s->state[k];
+    }
+
+    nb--;
+    s->visited = true;
+    for (const auto &successor : s->successors) {
+      if (!successor.first->visited) {
+        nb++;
+        TabBDDNodes(successor.first, nb);
+      }
+    }
+  }
+}
+
+void ModularSOG::PrintFullInformation(int nb_subnets) {
+  std::cout << "\n\nGRAPH SIZE : \n";
+  std::cout << "\n\tNB MARKING : " << nb_marking;
+  std::cout << "\n\tNB NODES : " << nb_states;
+  std::cout << "\n\tNB ARCS : " << nb_arcs << '\n';
+
+  std::cout << " \n\nCOMPLETE INFORMATION ?(y/n)\n";
+  char c = 0;
+  std::cin >> c;
+
+  tmp_table = new bdd[nb_states * nb_subnets];
+  size_t bdd_node = 1;
+  TabBDDNodes(initial_state, bdd_node);
+  std::cout << "NB BDD NODES : "
+            << bdd_anodecount(tmp_table, nb_states * nb_subnets) << '\n';
+  InitVisit(initial_state, 1);
+  if (c == 'y' || c == 'Y') {
+    size_t node = 1;
+    PrintGraph(initial_state, node);
+  }
+}
+
+void ModularSOG::PrintGraph(ModularAggregate *s, size_t &nb) const {
+  if (nb <= nb_states) {
+    std::cout << "\nSTATE NUMBER " << nb << " : \n";
+    s->visited = true;
+    PrintSuccessors(s);
+    getchar();
+    PrintPredecessors(s);
+    getchar();
+    for (const auto &successor : s->successors) {
+      if (!successor.first->visited) {
+        nb++;
+        PrintGraph(successor.first, nb);
+      }
+    }
+  }
+}
+
+void ModularSOG::Reset() {
+  current_state = nullptr;
+  nb_arcs = nb_marking = nb_states = 0;
+}
+
+void ModularSOG::PrintSuccessors(const ModularAggregate *s) const {
+  std::cout << *s << '\n';
+
+  if (s->boucle) {
+    std::cout << "\n\tON BOUCLE DESSUS AVEC EPSILON\n";
+  }
+
+  if (s->blocage) {
+    std::cout << "\n\tEXISTENCE D'UN ETAT BLOCANT\n";
+  }
+
+  std::cout << "\n\tSES SUCCESSEURS SONT  ( " << s->successors.size()
+            << " ) :\n\n";
+  getchar();
+
+  for (const auto &successor : s->successors) {
+    std::cout << " \t---" << successor.second << " -->";
+    std::cout << *(successor.first) << '\n';
+    getchar();
+  }
+}
+
+void ModularSOG::PrintPredecessors(const ModularAggregate *s) const {
+  std::cout << "\n\tSES PREDESCESSEURS SONT  ( " << s->predecessors.size()
+            << " ) :\n\n";
+  getchar();
+  for (const auto &predecessor : s->predecessors) {
+    std::cout << " \t---" << predecessor.second << " -->";
+    std::cout << *predecessor.first << '\n';
+    getchar();
+  }
+}
+
+void ModularSOG::IncrementNbArcs() {
+  nb_arcs++;
+}
diff --git a/src/modular_sog.hpp b/src/modular_sog.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..c10b1a06d8f729ca28c7dcf163f204062262e903
--- /dev/null
+++ b/src/modular_sog.hpp
@@ -0,0 +1,107 @@
+#ifndef MODULAR_SOG_H_
+#define MODULAR_SOG_H_
+
+#include <set>
+#include <vector>
+
+#include "modular_aggregate.hpp"
+
+// define types
+typedef std::set<int> Set;
+typedef std::vector<ModularAggregate *> ModularAggregates;
+
+class ModularSOG {
+ public:
+  ModularSOG() = default;
+
+  ModularAggregates GONodes;
+
+  ModularAggregate *initial_state{nullptr};
+  ModularAggregate *current_state{nullptr};
+  size_t nb_states{0};
+  size_t nb_marking{0};
+  size_t nb_arcs{0};
+
+  /**
+   * Reset the graph
+   */
+  void Reset();
+
+  /**
+   * Find a state
+   * @return
+   */
+  ModularAggregate *FindState(const ModularAggregate *s) const;
+
+  /**
+   * Find a modular state
+   * @return
+   */
+  ModularAggregate *FindState2(const ModularAggregate *s,
+                               Set states_involved) const;
+
+  /**
+   * Print successors of an state
+   */
+  void PrintSuccessors(const ModularAggregate *s) const;
+
+  /**
+   * Print predecessors of an state
+   */
+  void PrintPredecessors(const ModularAggregate *s) const;
+
+  /**
+   * Reset the visit flag of all states
+   * @param s
+   * @param nb
+   */
+  void InitVisit(ModularAggregate *s, size_t nb);
+
+  /**
+   * Fill the intermediate table to calculate the size of the graph
+   * @param s
+   * @param nb
+   */
+  void TabBDDNodes(ModularAggregate *s, size_t &nb);
+
+  /**
+   * Insert a state
+   * @param s
+   */
+  void InsertState(ModularAggregate *s);
+
+  /**
+   * Add an arc
+   * TODO: check the implementation of this
+   */
+  void IncrementNbArcs();
+
+  /**
+   * Add arcs to the graph
+   * @param pred
+   * @param succ
+   * @param trans
+   */
+  void AddArc(ModularAggregate *pred, ModularAggregate *succ,
+              const char *trans);
+
+  /**
+   * Print the graph
+   * @param s
+   * @param nb
+   */
+  void PrintGraph(ModularAggregate *s, size_t &nb) const;
+
+  /**
+   * Set the initial state of this graph
+   * @param s
+   */
+  void set_initial_state(ModularAggregate *s);
+
+  /**
+   * Print full information of the graph
+   * @param nb_subnets
+   */
+  void PrintFullInformation(int nb_subnets);
+};
+#endif  // MODULAR_SOG_H_
diff --git a/src/petri_net_sog.cpp b/src/petri_net_sog.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..3b8880d8b5f12864ae772c5ba0c6b357de4a91f3
--- /dev/null
+++ b/src/petri_net_sog.cpp
@@ -0,0 +1,840 @@
+#include "petri_net_sog.hpp"
+
+#include <iostream>
+#include <map>
+#include <set>
+#include <stack>
+#include <string>
+#include <utility>
+#include <vector>
+
+#include "Net.hpp"
+#include "bvec.h"
+
+using namespace std;
+
+// BDD initial values
+constexpr int BDD_INITIAL_NUM_NODES_ = 1000000;
+constexpr int BDD_SIZE_CACHES = 1000000;
+
+// TODO: why these global variables ?
+int nb_iter;
+int it_ext, it_int;
+int max_bdd_nodes;
+bdd *tab_meta;
+int nb_meta_state;
+double old_size;
+const vector<Place> *v_places = nullptr;
+
+/**
+ * Handler for errors in the BDD package
+ * @param err_code
+ */
+void BDDErrorHandler(int err_code) {
+  bdd_default_errhandler(err_code);
+}
+
+/**
+ * Handler to convert the FDD integer identifier into something readable by the
+ * end user
+ * @param o
+ * @param var
+ */
+void PrintHandler(ostream &o, int var) {
+  o << (*v_places)[var / 2].name;
+  if (var % 2) {
+    o << "_p";
+  }
+}
+
+bool operator<(const pair<bdd, bdd> &r, const pair<bdd, bdd> &s) {
+  return r.first.id() == s.first.id();
+}
+
+int choisir(Set cov, Set fire) {
+  auto it = fire.begin();
+  while (it != fire.end()) {
+    if (cov.find(*it) == cov.end()) {
+      return *it;
+    }
+    ++it;
+  }
+
+  return -1;
+}
+
+void reinit_cycle(const vector<int> &sw, map<int, int> &trans_obs) {
+  for (auto i : sw) {
+    if (trans_obs[i] > 0) {
+      trans_obs[i] = 2;
+    }
+  }
+}
+
+/*****************************************************************************/
+/*                                Class Trans                                */
+/*****************************************************************************/
+Trans::Trans(const bdd &v, bddPair *p, const bdd &r, const bdd &pr,
+             const bdd &pre, const bdd &post)
+    : var(v), pair(p), precond(pre), postcond(post), rel(r), prerel(pr) {}
+
+bdd Trans::operator()(const bdd &op) const {
+  bdd res = bdd_relprod(op, rel, var);
+  res = bdd_replace(res, pair);
+  return res;
+}
+
+bdd Trans::operator[](const bdd &op) const {
+  bdd res = bdd_relprod(op, prerel, var);
+  res = bdd_replace(res, pair);
+  return res;
+}
+
+/*****************************************************************/
+/*                         Class RdPBDD                          */
+/*****************************************************************/
+
+PetriNetSOG::PetriNetSOG(const net &petri_net, const map<int, int> &obs_trans,
+                         Set non_obs_trans, const int bound, const bool init)
+    : nb_places(petri_net.places.size()) {
+  bvec *v = new bvec[nb_places];
+  bvec *vp = new bvec[nb_places];
+  bvec *prec = new bvec[nb_places];
+  bvec *postc = new bvec[nb_places];
+
+  int *idv = new int[nb_places];
+  int *idvp = new int[nb_places];
+  int *nbbit = new int[nb_places];
+
+  if (!init) {
+    bdd_init(BDD_INITIAL_NUM_NODES_, BDD_SIZE_CACHES);
+  }
+
+  // the error handler
+  bdd_error_hook(BDDErrorHandler);
+
+  // add petri net transitions to the set of observable transitions
+  transitions = petri_net.transitions;
+  for (auto t : obs_trans) {
+    observables.insert(t.first);
+  };
+
+  // since non_observables is passed by value, we can move it
+  non_observables = std::move(non_obs_trans);
+
+  transitions_names = petri_net.transitionName;
+
+  int domain = 0;
+  for (const auto &place : petri_net.places) {
+    if (place.hasCapacity()) {
+      domain = place.capacity + 1;
+    } else {
+      domain = bound + 1;  // the default domain
+    }
+
+    // variables are created one by one (implying contiguous binary variables)
+    fdd_extdomain(&domain, 1);
+    fdd_extdomain(&domain, 1);
+  }
+
+  // bvec
+  current_var = bdd_true();
+  for (int i = 0; i < nb_places; i++) {
+    nbbit[i] = fdd_varnum(2 * i);
+    v[i] = bvec_varfdd(2 * i);
+    vp[i] = bvec_varfdd((2 * i) + 1);
+    current_var = current_var & fdd_ithset(2 * i);
+  }
+
+  // initial marking
+  m0 = bdd_true();
+  int offset = 0;
+  for (const auto &place : petri_net.places) {
+    m0 = m0 & fdd_ithvar(2 * offset, place.marking);
+    offset++;
+  }
+
+  // place names
+  v_places = &petri_net.places;
+  fdd_strm_hook(PrintHandler);
+
+  /* Transition relation */
+  for (auto t = petri_net.transitions.begin(); t != petri_net.transitions.end();
+       ++t) {
+    bdd rel = bdd_true();
+    bdd var = bdd_true();
+    bdd prerel = bdd_true();
+    bdd precond = bdd_true();
+    bdd postcond = bdd_true();
+
+    bddPair *p = bdd_newpair();
+
+    for (int i = 0; i < nb_places; i++) {
+      prec[i] = bvec_con(nbbit[i], 0);
+      postc[i] = bvec_con(nbbit[i], 0);
+    }
+
+    // arcs pre
+    set<int> adjacent_places;
+    for (auto it = t->pre.begin(); it != t->pre.end(); ++it) {
+      adjacent_places.insert(it->first);
+      prec[it->first] =
+          prec[it->first] + bvec_con(nbbit[it->first], it->second);
+    }
+
+    // arcs post
+    for (auto it = t->post.begin(); it != t->post.end(); ++it) {
+      adjacent_places.insert(it->first);
+      postc[it->first] =
+          postc[it->first] + bvec_con(nbbit[it->first], it->second);
+    }
+
+    int np = 0;
+    for (const auto adjacent_place : adjacent_places) {
+      idv[np] = 2 * adjacent_place;
+      idvp[np] = 2 * adjacent_place + 1;
+      var = var & fdd_ithset(2 * adjacent_place);
+
+      // image
+      // precondition
+      rel = rel & (v[adjacent_place] >= prec[adjacent_place]);
+      precond = precond & (v[adjacent_place] >= prec[adjacent_place]);
+      // postcondition
+      rel = rel &
+            (vp[adjacent_place] == (v[adjacent_place] - prec[adjacent_place] +
+                                    postc[adjacent_place]));
+
+      // pre-image
+      // precondition
+      prerel = prerel & (v[adjacent_place] >= postc[adjacent_place]);
+      // postcondition
+      postcond = postcond & (v[adjacent_place] >= postc[adjacent_place]);
+      prerel =
+          prerel &
+          (vp[adjacent_place] ==
+           (v[adjacent_place] - postc[adjacent_place] + prec[adjacent_place]));
+
+      // capacity
+      if (petri_net.places[adjacent_place].hasCapacity()) {
+        rel = rel & (vp[adjacent_place] <=
+                     bvec_con(nbbit[adjacent_place],
+                              petri_net.places[adjacent_place].capacity));
+      }
+      np++;
+    }
+
+    fdd_setpairs(p, idvp, idv, np);
+    relation.emplace_back(var, p, rel, prerel, precond, postcond);
+  }
+
+  // remove vectors
+  delete[] v;
+  delete[] vp;
+  delete[] prec;
+  delete[] postc;
+  delete[] idv;
+  delete[] idvp;
+  delete[] nbbit;
+}
+
+bdd PetriNetSOG::ReachableBDD1() const {
+  bdd m2 = m0;
+  nb_iter = 0;
+  max_bdd_nodes = bdd_nodecount(m0);
+
+  bdd m1;
+  while (m1 != m2) {
+    m1 = m2;
+    for (const auto &trans : relation) {
+      m2 = trans(m2) | m2;
+    }
+
+    const int current_size = bdd_nodecount(m2);
+    if (max_bdd_nodes < current_size) {
+      max_bdd_nodes = current_size;
+    }
+
+    nb_iter++;
+  }
+
+  return m2;
+}
+
+bdd PetriNetSOG::ReachableBDD2() const {
+  bdd m2 = m0;
+
+  nb_iter = 0;
+  max_bdd_nodes = bdd_nodecount(m0);
+  bdd m1;
+  while (m1 != m2) {
+    m1 = m2;
+
+    bdd reached;
+    for (const auto &trans : relation) {
+      reached = trans(m2) | reached;
+    }
+    m2 = m2 | reached;
+
+    const int current_size = bdd_nodecount(m2);
+    if (max_bdd_nodes < current_size) {
+      max_bdd_nodes = current_size;
+    }
+
+    nb_iter++;
+  }
+
+  return m2;
+}
+
+bdd PetriNetSOG::ReachableBDD3() const {
+  bdd new_bdd;
+  bdd from = m0;
+  bdd reached = m0;
+  nb_iter = 0;
+
+  do {
+    nb_iter++;
+
+    bdd succ;
+    for (const auto &trans : relation) {
+      succ = trans(from) | succ;
+    }
+
+    new_bdd = succ - reached;
+    from = new_bdd;
+    reached = reached | new_bdd;
+  } while (new_bdd != bddfalse);
+
+  return reached;
+}
+
+bdd PetriNetSOG::AccessibleEpsilon(const bdd &from_s) const {
+  bdd m1;
+  bdd m2 = from_s;
+
+  do {
+    m1 = m2;
+    for (const int i : non_observables) {
+      m2 = relation[i](m2) | m2;
+    }
+  } while (m1 != m2);
+
+  return m2;
+}
+
+bdd PetriNetSOG::AccessibleEpsilon2(const bdd &from_s) const {
+  bdd new_bdd;
+  bdd from = from_s;
+  bdd reached = from_s;
+
+  do {
+    bdd succ;
+    for (const int i : non_observables) {
+      succ = relation[i](from) | succ;
+    }
+
+    new_bdd = succ - reached;
+    from = new_bdd;
+    reached = reached | new_bdd;
+  } while (new_bdd != bddfalse);
+
+  return reached;
+}
+
+bdd PetriNetSOG::StepForward(const bdd &from) const {
+  bdd res = from;
+  for (const int i : non_observables) {
+    res = res | relation[i](res);
+  }
+  return res;
+}
+
+bdd PetriNetSOG::StepForward2(const bdd &from) const {
+  bdd res;
+  for (const int i : non_observables) {
+    res = res | relation[i](from);
+  }
+  return res;
+}
+
+bdd PetriNetSOG::StepBackward(const bdd &from) const {
+  bdd res = from;
+
+  for (const auto &t : relation) {
+    res = res | t[res];
+  }
+
+  return res;
+}
+
+pair<int, bdd> PetriNetSOG::StepBackward1(const bdd &from,
+                                          const Aggregate *aggr) const {
+  pair<int, bdd> res;
+
+  for (const auto t : non_observables) {
+    bdd succ = (relation[t])[from];
+
+    // function that returns the preceding bdd with the transition t
+    if ((succ != bdd_false()) & ((succ &= aggr->state) != bdd_false())) {
+      res.first = t;
+      res.second = succ;
+      break;
+    }
+  }
+
+  return res;
+}
+
+bdd PetriNetSOG::StepBackward2(const bdd &from) const {
+  bdd res;
+
+  for (const auto &t : relation) {
+    res = res | t[from];
+  }
+
+  return res;
+}
+
+bdd PetriNetSOG::GetSuccessor(const bdd &from, const int t) const {
+  return relation[t](from);
+}
+
+Set PetriNetSOG::FirableObs(const bdd &state) const {
+  Set res;
+
+  for (int i : observables) {
+    if (relation[i](state) != bddfalse) {
+      res.insert(i);
+    }
+  }
+
+  return res;
+}
+
+set<chem> PetriNetSOG::ObservablePaths(SOG &g, map<int, int> trans_obs) const {
+  set<chem> cto;
+  vector<int> sw;
+  Set cov;
+  Stack st;
+  Set fire;
+
+  tab_meta = new bdd[100000];
+  nb_meta_state = 0;
+  max_bdd_nodes = 0;
+  nb_iter = 0;
+  it_ext = it_int = 0;
+
+  bool ancien = true;
+
+  // construction of the first aggregate
+  auto *c = new Aggregate;
+  {
+    bdd complete_meta_state = AccessibleEpsilon(m0);
+    fire = FirableObs(complete_meta_state);
+    c->state = complete_meta_state;
+    tab_meta[nb_meta_state] = c->state;
+    nb_meta_state++;
+    old_size = bdd_nodecount(c->state);
+    st.emplace(Couple(c, complete_meta_state), fire);
+  }
+
+  g.set_initial_state(c);
+  g.InsertState(c);
+  g.nb_marking += bdd_pathcount(c->state);
+
+  while (!st.empty()) {
+    StackElt e = st.top();
+    st.pop();
+    nb_meta_state--;
+
+    // if there are firable transitions
+    if (!e.second.empty()) {
+      // choose a transition from the firable transitions
+      int t = choisir(cov, e.second);
+      if (t != -1) {
+        ancien = false;
+        trans_obs[t]--;
+        sw.push_back(t);
+        if (trans_obs[t] == 0) {
+          cov.insert(t);
+          if (cov.size() == observables.size()) {
+            cto.insert(sw);
+            return cto;
+          }
+        }
+      } else {
+        t = *e.second.begin();
+        sw.push_back(t);
+      }
+
+      e.second.erase(t);
+      st.push(e);
+
+      // double nbnode;
+      {
+        auto *reached_class = new Aggregate;
+        bdd complete_meta_state =
+            AccessibleEpsilon(GetSuccessor(e.first.second, t));
+        reached_class->state = complete_meta_state;
+        Aggregate *pos = g.FindState(reached_class);
+
+        // new aggregate
+        if (!pos) {
+          fire = FirableObs(complete_meta_state);
+
+          nb_meta_state++;
+
+          // add transitions
+          e.first.first->successors.insert(e.first.first->successors.begin(),
+                                           Edge(reached_class, t));
+          reached_class->predecessors.insert(
+              reached_class->predecessors.begin(), Edge(e.first.first, t));
+          g.IncrementNbArcs();
+
+          // add the aggregate
+          g.InsertState(reached_class);
+          if (fire.empty()) {
+            if (!ancien) {
+              cto.insert(sw);
+            }
+            reinit_cycle(sw, trans_obs);
+
+            sw.pop_back();
+          } else {
+            st.emplace(Couple(reached_class, complete_meta_state), fire);
+          }
+        } else {  // aggregate already exists
+          if (!ancien) {
+            cto.insert(sw);
+          }
+          reinit_cycle(sw, trans_obs);
+
+          sw.pop_back();
+          ancien = true;
+          delete reached_class;
+          e.first.first->successors.insert(e.first.first->successors.begin(),
+                                           Edge(pos, t));
+          pos->predecessors.insert(pos->predecessors.begin(),
+                                   Edge(e.first.first, t));
+          g.IncrementNbArcs();
+        }
+      }
+    } else {
+      sw.pop_back();
+      ancien = true;
+    }
+  }
+
+  return cto;
+}
+
+stack<Couple> PetriNetSOG::RecherchePointsEntree(chem path,
+                                                 const SOG &g) const {
+  Couple p;
+  stack<Couple> pt_entr;
+
+  bdd entree = m0;
+  Aggregate *agr = g.initial_state;
+
+  p.first = agr;
+  p.second = entree;
+  pt_entr.push(p);
+
+  for (auto k = path.begin(); k != path.end() - 1; ++k) {
+    const int t = *k;
+    entree = relation[t](p.first->state);
+    p.second = entree;
+
+    for (const auto &succ : agr->successors) {
+      if (succ.second == t) {
+        agr = succ.first;
+        break;
+      }
+    }
+
+    p.first = agr;
+    pt_entr.push(p);
+  }
+
+  return pt_entr;
+}
+
+chem PetriNetSOG::AbstractPaths(chem path, const SOG &g) const {
+  bdd source;
+  vector<int> abstract_paths;
+  stack<Couple> point_entree = RecherchePointsEntree(path, g);
+
+  int i = 0;
+  while (!point_entree.empty()) {
+    int trans = *(path.end() - 1);
+    const Couple agr_entree = point_entree.top();
+    point_entree.pop();
+
+    Aggregate *agr = agr_entree.first;
+    const bdd target = agr_entree.second;
+
+    if (i == 0) {
+      source = FrontiereNodes1(agr->state, trans);
+    }
+
+    if (i == 1) {
+      source = (relation[trans])[source];
+    }
+
+    vector<int> path_aggregate = SubPathAggregate(&source, target, agr);
+    abstract_paths.insert(abstract_paths.begin(), trans);
+    abstract_paths.insert(abstract_paths.begin(), path_aggregate.begin(),
+                          path_aggregate.end());
+    path.pop_back();
+    i = 1;
+  }
+
+  return abstract_paths;
+}
+
+vector<int> PetriNetSOG::SubPathAggregate(bdd *source, const bdd &target,
+                                          const Aggregate *aggr) const {
+  vector<int> path;
+
+  bdd current_state = *source;
+  while ((target & current_state) == bdd_false()) {
+    pair<int, bdd> couple = StepBackward1(current_state, aggr);
+    path.insert(path.begin(), couple.first);
+    current_state = couple.second;
+  }
+
+  *source = current_state;
+  return path;
+}
+
+bdd PetriNetSOG::CanonizeR(bdd s, unsigned int i) {
+  bdd s1, s2;
+  do {
+    it_ext++;
+    s1 = s - bdd_nithvar(2 * i);
+    s2 = s - s1;
+
+    // TODO: extract the loop into a function
+    if (s1 != bddfalse && s2 != bddfalse) {
+      bdd front = s1;
+      bdd reached = s1;
+      do {
+        it_int++;
+        front = StepForward(front) - reached;
+        reached = reached | front;
+        s2 = s2 - front;
+      } while (front != bddfalse && s2 != bddfalse);
+    }
+
+    if (s1 != bddfalse && s2 != bddfalse) {
+      bdd front = s2;
+      bdd reached = s2;
+      do {
+        it_int++;
+        front = StepForward(front) - reached;
+        reached = reached | front;
+        s1 = s1 - front;
+      } while (front != bddfalse && s1 != bddfalse);
+    }
+
+    s = s1 | s2;
+    i++;
+  } while (i < nb_places && (s1 == bddfalse || s2 == bddfalse));
+
+  if (i >= nb_places) {
+    return s;
+  }
+
+  return (CanonizeR(s1, i) | CanonizeR(s2, i));
+}
+
+bool PetriNetSOG::SetBlocking(const bdd &s) const {
+  bdd mt = bddtrue;
+
+  for (const auto &i : relation) {
+    mt = mt & !(i.precond);
+  }
+
+  return ((s & mt) != bddfalse);
+}
+
+// TODO: check again this
+bool PetriNetSOG::SetDiv(const bdd &s) const {
+  bdd reached = s;
+
+  do {
+    bdd from = reached;
+    for (const int trans : non_observables) {
+      bdd to = relation[trans](reached);
+      reached = reached | to;  // Reached=To ???
+    }
+
+    if (reached == from) {
+      return true;
+    }
+  } while (reached != bddfalse);
+
+  return false;
+}
+
+bdd PetriNetSOG::FrontiereNodes(const bdd &from) const {
+  bdd res = bddfalse;
+  for (const int observable : observables) {
+    res = res | (from & relation[observable].precond);
+  }
+
+  return res;
+}
+
+bdd PetriNetSOG::FrontiereNodes1(const bdd &from, const int t) const {
+  return bddfalse | (from & relation[t].precond);
+}
+
+void PetriNetSOG::GeneralizedSyncProduct(ModularSOG &gv, const int nb_subnets,
+                                         PetriNetSOG *subnets[]) {
+  // TODO: what is that?
+  int pos_trans(TRANSITIONS, string);
+
+  tab_meta = new bdd[1000000];
+  nb_meta_state = 0;
+  max_bdd_nodes = 0;
+  nb_meta_state = 0;
+
+  auto complete_meta_state = new bdd[nb_subnets];
+  auto fire = new Set[nb_subnets];
+  auto meta_state = new ModularAggregate;
+
+  for (int k = 0; k < nb_subnets; k++) {
+    complete_meta_state[k] = subnets[k]->AccessibleEpsilon(subnets[k]->m0);
+    fire[k] = subnets[k]->FirableObs(complete_meta_state[k]);
+    meta_state->state.push_back(
+        subnets[k]->CanonizeR(complete_meta_state[k], 0));
+    tab_meta[nb_meta_state] = meta_state->state[k];
+    nb_meta_state++;
+  }
+
+  old_size = bdd_anodecount(tab_meta, nb_meta_state);
+  meta_state->blocage = true;
+  for (int k = 0; ((k < nb_subnets) && (meta_state->blocage)); k++) {
+    meta_state->blocage =
+        meta_state->blocage && subnets[k]->SetBlocking(complete_meta_state[k]);
+  }
+  meta_state->boucle = false;
+
+  for (int k = 0; ((k < nb_subnets) && (!meta_state->boucle)); k++) {
+    meta_state->boucle =
+        meta_state->boucle || subnets[k]->SetDiv(complete_meta_state[k]);
+  }
+
+  gv.set_initial_state(meta_state);
+  gv.InsertState(meta_state);
+  nb_meta_state++;
+
+  ModularStack st;
+  st.emplace(ModularCouple(meta_state, complete_meta_state), fire);
+  do {
+    nb_iter++;
+    StackModularElt e = st.top();
+    st.pop();
+    for (int k = 0; k < nb_subnets; k++) {
+      while (!e.second[k].empty()) {
+        int t = *e.second[k].begin();
+        e.second[k].erase(t);
+        bool ok = true;
+        Set concerned_subnets;
+        concerned_subnets.insert(k);
+        string tmp = subnets[k]->transitions[t].name;
+        for (int j = 0; j < nb_subnets; j++) {
+          if (j != k) {
+            int num = subnets[j]->transitions_names[tmp];
+            const int pos = pos_trans(subnets[j]->transitions, tmp);
+            if (pos != -1 && subnets[j]->observables.find(num) !=
+                                 subnets[j]->observables.end()) {
+              concerned_subnets.insert(j);
+              if (e.second[j].find(num) == e.second[j].end()) {
+                ok = false;
+              } else
+                e.second[j].erase(num);
+            }
+          }
+        }
+        if (ok) {
+          complete_meta_state = new bdd[nb_subnets];
+          fire = new Set[nb_subnets];
+          meta_state = new ModularAggregate;
+
+          for (int j = 0; j < nb_subnets; j++) {
+            if (concerned_subnets.find(j) == concerned_subnets.end()) {
+              complete_meta_state[j] = e.first.second[j];
+              meta_state->state.push_back(e.first.first->state[j]);
+            } else {
+              complete_meta_state[j] =
+                  subnets[j]->AccessibleEpsilon(subnets[j]->GetSuccessor(
+                      e.first.second[j], subnets[j]->transitions_names[tmp]));
+              meta_state->state.push_back(
+                  subnets[j]->CanonizeR(complete_meta_state[j], 0));
+              tab_meta[nb_meta_state] = meta_state->state[k];
+              nb_meta_state++;
+            }
+
+            fire[j] = subnets[j]->FirableObs(complete_meta_state[j]);
+          }
+
+          ModularAggregate *pos = gv.FindState(meta_state);
+          if (!pos) {
+            old_size = bdd_anodecount(tab_meta, nb_meta_state);
+            // Computes the deadlock and loop attributes
+            meta_state->blocage = true;
+            for (int j = 0; ((j < nb_subnets) && (meta_state->blocage)); j++) {
+              meta_state->blocage &=
+                  subnets[j]->SetBlocking(complete_meta_state[j]);
+            }
+
+            meta_state->boucle = false;
+            for (int j = 0; ((j < nb_subnets) && (!meta_state->boucle)); j++) {
+              meta_state->boucle |= subnets[j]->SetDiv(complete_meta_state[j]);
+            }
+
+            st.emplace(ModularCouple(meta_state, complete_meta_state), fire);
+            e.first.first->successors.insert(e.first.first->successors.begin(),
+                                             ModularEdge(meta_state, tmp));
+            meta_state->predecessors.insert(meta_state->predecessors.begin(),
+                                            ModularEdge(e.first.first, tmp));
+            gv.IncrementNbArcs();
+            gv.InsertState(meta_state);
+          } else {
+            e.first.first->successors.insert(e.first.first->successors.begin(),
+                                             ModularEdge(pos, tmp));
+            pos->predecessors.insert(pos->predecessors.begin(),
+                                     ModularEdge(e.first.first, tmp));
+            gv.IncrementNbArcs();
+            delete meta_state;
+          }
+        }
+      }
+    }
+  } while (!st.empty());
+}
+
+bdd PetriNetSOG::EmersonLey(const bdd &state, const bool trace) const {
+  bdd b = state;
+  bdd old_b;
+  bdd old_d;
+
+  const bdd fair = bdd_ithvar(2 * nb_places - 1);
+  do {
+    old_b = b;
+    bdd d = b & fair;
+    do {
+      old_d = d;
+      bdd inter = b & StepForward2(d);
+      d = d | inter;
+    } while (old_d != d);
+
+    if (trace) {
+      b = b & StepBackward2(d);
+    }
+  } while (b != old_b);
+
+  return b;
+}
diff --git a/src/petri_net_sog.hpp b/src/petri_net_sog.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..1b3e0aec51dab97541fdf4f7b109cf2158cd3064
--- /dev/null
+++ b/src/petri_net_sog.hpp
@@ -0,0 +1,131 @@
+#ifndef PETRI_NET_SOG_H_
+#define PETRI_NET_SOG_H_
+#include <stack>
+#include <vector>
+
+#include "Net.hpp"
+#include "bdd.h"
+#include "modular_aggregate.hpp"
+#include "modular_sog.hpp"
+#include "sog.hpp"
+
+// type definitions
+typedef std::vector<int> chem;
+
+// structures useful for the synchronised product of n observation graphs.
+typedef pair<Aggregate*, bdd> Couple;
+typedef std::pair<ModularAggregate*, bdd*> ModularCouple;
+
+typedef pair<Couple, Set> StackElt;
+typedef std::pair<ModularCouple, Set*> StackModularElt;
+
+typedef stack<StackElt> Stack;
+typedef std::stack<StackModularElt> ModularStack;
+
+class Trans {
+ public:
+  Trans(const bdd& var, bddPair* pair, const bdd& rel, const bdd& prerel,
+        const bdd& precond, const bdd& postcond);
+
+  /**
+   * Forward firing
+   * @param op
+   * @return
+   */
+  bdd operator()(const bdd& op) const;
+
+  /**
+   * Backward firing
+   * @param op
+   * @return
+   */
+  bdd operator[](const bdd& op) const;
+
+  bdd var;
+  bddPair* pair;
+  bdd precond;
+  bdd postcond;
+  bdd rel;
+  bdd prerel;
+};
+
+class PetriNetSOG {
+ private:
+  std::vector<Transition> transitions;
+  Set observables;
+  Set non_observables;
+  std::map<std::string, int> transitions_names;
+  size_t nb_places;
+
+ public:
+  // initial marking
+  bdd m0;
+  bdd current_var;
+  std::vector<Trans> relation;
+
+  PetriNetSOG(const net&, const map<int, int>& obs_trans, Set non_obs_trans,
+              int bound = 32, bool init = false);
+  ~PetriNetSOG() = default;
+
+  // reachability functions
+  bdd ReachableBDD1() const;
+  bdd ReachableBDD2() const;
+  bdd ReachableBDD3() const;
+
+  // transitive closure on unobserved transitions
+  bdd AccessibleEpsilon(const bdd& from_s) const;
+  bdd AccessibleEpsilon2(const bdd& from_s) const;
+
+  // step forward using non-observable transitions
+  bdd StepForward(const bdd& from) const;
+  bdd StepForward2(const bdd& from) const;
+
+  // step backward
+  bdd StepBackward(const bdd& from) const;
+  std::pair<int, bdd> StepBackward1(const bdd& from,
+                                    const Aggregate* aggr) const;
+  bdd StepBackward2(const bdd& from) const;
+
+  /**
+   * Returns the successor of a state from a transition
+   * @param from
+   * @param t
+   * @return
+   */
+  bdd GetSuccessor(const bdd& from, int t) const;
+
+  /**
+   * Return the set of observable transitions that are firable from a state
+   * @param state
+   * @return
+   */
+  Set FirableObs(const bdd& state) const;
+
+  // Extract observable paths
+  std::set<chem> ObservablePaths(SOG& g, std::map<int, int> trans_obs) const;
+
+  // SOG generation
+  std::stack<Couple> RecherchePointsEntree(chem path, const SOG& g) const;
+
+  // Generate abstract paths
+  chem AbstractPaths(chem path, const SOG& g) const;
+
+  // Find a path from a source state to a target state in an aggregate
+  std::vector<int> SubPathAggregate(bdd* source, const bdd& target,
+                                    const Aggregate* aggr) const;
+
+  bdd CanonizeR(bdd s, unsigned int i);
+
+  bool SetBlocking(const bdd& s) const;
+  bool SetDiv(const bdd& s) const;
+
+  bdd FrontiereNodes(const bdd& from) const;
+  bdd FrontiereNodes1(const bdd& from, int t) const;
+
+  void GeneralizedSyncProduct(ModularSOG& gv, int nb_subnets,
+                              PetriNetSOG* subnets[]);
+
+  bdd EmersonLey(const bdd& state, bool trace) const;
+};
+
+#endif  // PETRI_NET_SOG_H_
diff --git a/src/sog.cpp b/src/sog.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f13e8d200a1f8f3d2c54b4349e0a7553d107b0f7
--- /dev/null
+++ b/src/sog.cpp
@@ -0,0 +1,151 @@
+
+#include "sog.hpp"
+
+#include <fstream>
+
+bdd *tab;
+
+void SOG::set_initial_state(Aggregate *s) {
+  current_state = initial_state = s;
+}
+
+Aggregate *SOG::FindState(const Aggregate *s) const {
+  for (auto *GONode : GONodes) {
+    if (s->state.id() == GONode->state.id()) {
+      return GONode;
+    }
+  }
+  return nullptr;
+}
+
+void SOG::InsertState(Aggregate *s) {
+  s->visited = false;
+  this->GONodes.push_back(s);
+  nb_states++;
+}
+
+int SOG::NbBddNode(Aggregate *s, size_t &nb) {
+  if (!s->visited) {
+    tab[nb - 1] = s->state;
+    s->visited = true;
+    const int bdd_node = bdd_nodecount(s->state);
+    int size_succ = 0;
+    for (const auto &successor : s->successors) {
+      if (!successor.first->visited) {
+        nb++;
+        size_succ += NbBddNode(successor.first, nb);
+      }
+    }
+    return size_succ + bdd_node;
+  }
+
+  return 0;
+}
+
+void SOG::PrintCompleteInformation() {
+  std::cout << "\n\nGRAPH SIZE :";
+  std::cout << "\n\tNB MARKING : " << nb_marking;
+  std::cout << "\n\tNB NODES : " << nb_states;
+  std::cout << "\n\tNB ARCS : " << nb_arcs << "\n";
+
+  tab = new bdd[static_cast<int>(nb_states)];
+
+  size_t node = 1;
+  NbBddNode(initial_state, node);
+  std::cout << "\tNB BDD NODES : "
+            << bdd_anodecount(tab, static_cast<int>(nb_states)) << "\n";
+  InitVisit(initial_state, 1);
+}
+
+void SOG::InitVisit(Aggregate *s, size_t nb) {
+  if (nb <= nb_states) {
+    s->visited = false;
+    for (const auto &successor : s->successors) {
+      if (successor.first->visited) {
+        nb++;
+        InitVisit(successor.first, nb);
+      }
+    }
+  }
+}
+
+void SOG::PrintGraph(Aggregate *s, size_t &nb) {
+  if (nb <= nb_states) {
+    std::cout << "\nSTATE NUMBER " << nb << " : \n";
+    s->visited = true;
+    PrintSuccessors(s);
+    getchar();
+    PrintPredecessors(s);
+    getchar();
+
+    for (const auto &successor : s->successors) {
+      if (!successor.first->visited) {
+        nb++;
+        PrintGraph(successor.first, nb);
+      }
+    }
+  }
+}
+
+void SOG::GenerateReachabilityGraphDotfile(const std::string &filename) const {
+  // Open file for writing
+  std::ofstream file("./result/" + filename + ".dot");
+
+  // Check if the file opened successfully
+  if (!file.is_open()) {
+    std::cerr << "Error: Could not open file for writing.\n";
+    return;
+  }
+
+  std::cout << "le graph contient : " << GONodes.size() << '\n';
+  file << "digraph reachab_graph {" << '\n';
+
+  for (const auto *GONode : GONodes) {
+    for (const auto &edge : GONode->successors) {
+      file << "ag_" << GONode->state.id() << " -> "
+           << "ag_" << edge.first->state.id() << " [ label =  \"t"
+           << edge.second + 1 << "\"]" << '\n';
+    }
+  }
+
+  file << "ag_" << initial_state->state.id() << "[initialstate]\n";
+  file << " }\n";
+  file.close();
+}
+
+void SOG::PrintSuccessors(const Aggregate *s) const {
+  std::cout << bddtable << s->state << '\n';
+  if (s->boucle) {
+    std::cout << "\n\tON BOUCLE DESSUS AVEC EPSILON\n";
+  }
+
+  if (s->blocage) {
+    std::cout << "\n\tEXISTENCE D'UN ETAT BLOCANT\n";
+  }
+
+  std::cout << "\n\tSES SUCCESSEURS SONT  ( " << s->successors.size()
+            << " ) :\n\n";
+  getchar();
+
+  for (const auto &successor : s->successors) {
+    std::cout << " \t- t" << successor.second << " ->" << bddtable
+              << successor.first->state << '\n';
+    getchar();
+  }
+}
+
+void SOG::PrintPredecessors(const Aggregate *s) const {
+  std::cout << "\n\tSES PREDESCESSEURS SONT  ( " << s->predecessors.size()
+            << " ) :\n\n";
+  getchar();
+
+  for (const auto &predecessor : s->predecessors) {
+    std::cout << " \t- t" << predecessor.second << " ->" << bddtable
+              << predecessor.first->state << '\n';
+    getchar();
+  }
+}
+
+void SOG::IncrementNbArcs() {
+  nb_arcs++;
+}
diff --git a/src/sog.hpp b/src/sog.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..094c4e5fcd0d7f094efc9a7d09f4097e49edcb25
--- /dev/null
+++ b/src/sog.hpp
@@ -0,0 +1,93 @@
+#ifndef SOG_H_
+#define SOG_H_
+
+#include <vector>
+
+#include "aggregate.hpp"
+
+typedef std::vector<Aggregate *> Aggregates;
+
+class SOG {
+ public:
+  SOG() = default;
+
+  Aggregates GONodes;
+
+  Aggregate *initial_state{nullptr};
+  Aggregate *current_state{nullptr};
+  size_t nb_states{0};
+  double nb_marking{0};
+  size_t nb_arcs{0};
+
+  /**
+   * Find a state
+   * @param s
+   * @return
+   */
+  Aggregate *FindState(const Aggregate *s) const;
+
+  /**
+   * Initialize the visit flag of all states
+   * @param s
+   * @param nb
+   */
+  void InitVisit(Aggregate *s, size_t nb);
+
+  /**
+   * Returns the number of BDD nodes of an state
+   * @param s
+   * @param nb
+   * @return
+   */
+  int NbBddNode(Aggregate *s, size_t &nb);
+
+  /**
+   * Print successors of an state
+   * @param s
+   */
+  void PrintSuccessors(const Aggregate *s) const;
+
+  /**
+   * Print predecessors of an state
+   * @param s
+   */
+  void PrintPredecessors(const Aggregate *s) const;
+
+  /**
+   * Add an arc
+   * TODO: please check the implementation of this function
+   */
+  void IncrementNbArcs();
+
+  /**
+   * Insert state
+   * @param s
+   */
+  void InsertState(Aggregate *s);
+
+  /**
+   * Set the initial state of this graph
+   * @param s
+   */
+  void set_initial_state(Aggregate *s);
+
+  /**
+   * Print the information of the graph
+   */
+  void PrintCompleteInformation();
+
+  /**
+   * Print the graph
+   * @param s
+   * @param nb
+   */
+  void PrintGraph(Aggregate *s, size_t &nb);
+
+  /**
+   * Print reachability graph using graphviz notation
+   * @param filename
+   */
+  void GenerateReachabilityGraphDotfile(const std::string &filename) const;
+};
+
+#endif  // SOG_H_