diff --git a/abclib.dsp b/abclib.dsp
index 4e027c495..db504b712 100644
--- a/abclib.dsp
+++ b/abclib.dsp
@@ -5211,6 +5211,10 @@ SOURCE=.\src\aig\gia\giaRex.c
 # End Source File
 # Begin Source File
 
+SOURCE=.\src\aig\gia\giaRrr.cpp
+# End Source File
+# Begin Source File
+
 SOURCE=.\src\aig\gia\giaSat3.c
 # End Source File
 # Begin Source File
diff --git a/src/aig/gia/gia.h b/src/aig/gia/gia.h
index 015a1da22..09e892160 100644
--- a/src/aig/gia/gia.h
+++ b/src/aig/gia/gia.h
@@ -1806,6 +1806,9 @@ extern Gia_Man_t *         Gia_ManTtoptCare( Gia_Man_t * p, int nIns, int nOuts,
 extern Gia_Man_t *         Gia_ManTransductionBdd( Gia_Man_t * pGia, int nType, int fMspf, int nRandom, int nSortType, int nPiShuffle, int nParameter, int fLevel, Gia_Man_t * pExdc, int fNewLine, int nVerbose );
 extern Gia_Man_t *         Gia_ManTransductionTt( Gia_Man_t * pGia, int nType, int fMspf, int nRandom, int nSortType, int nPiShuffle, int nParameter, int fLevel, Gia_Man_t * pExdc, int fNewLine, int nVerbose );
 
+/*=== giaRrr.cpp ===========================================================*/
+extern Gia_Man_t *         Gia_ManRrr( Gia_Man_t *pGia, int iSeed, int nWords, int nTimeout, int nSchedulerVerbose, int nOptimizerVerbose, int nAnalyzerVerbose, int nSimulatorVerbose, int nSatSolverVerbose, int fUseBddCspf, int fUseBddMspf, int nConflictLimit, int nSortType, int nOptimizerFlow, int nSchedulerFlow );
+
 /*=== giaCTas.c ===========================================================*/
 typedef struct Tas_Man_t_  Tas_Man_t;
 extern Tas_Man_t *         Tas_ManAlloc( Gia_Man_t * pAig, int nBTLimit );
diff --git a/src/aig/gia/giaRrr.cpp b/src/aig/gia/giaRrr.cpp
new file mode 100644
index 000000000..21e343337
--- /dev/null
+++ b/src/aig/gia/giaRrr.cpp
@@ -0,0 +1,31 @@
+#include "aig/gia/gia.h"
+
+#include "opt/rrr/rrr.h"
+#include "opt/rrr/rrrAbc.h"
+
+ABC_NAMESPACE_IMPL_START
+
+Gia_Man_t *Gia_ManRrr(Gia_Man_t *pGia, int iSeed, int nWords, int nTimeout, int nSchedulerVerbose, int nOptimizerVerbose, int nAnalyzerVerbose, int nSimulatorVerbose, int nSatSolverVerbose, int fUseBddCspf, int fUseBddMspf, int nConflictLimit, int nSortType, int nOptimizerFlow, int nSchedulerFlow) {
+  rrr::AndNetwork ntk;
+  ntk.Read<Gia_Man_t>(pGia, rrr::GiaReader<rrr::AndNetwork>);
+  rrr::Parameter Par;
+  Par.iSeed = iSeed;
+  Par.nWords = nWords;
+  Par.nTimeout = nTimeout;
+  Par.nSchedulerVerbose = nSchedulerVerbose;
+  Par.nOptimizerVerbose = nOptimizerVerbose;
+  Par.nAnalyzerVerbose = nAnalyzerVerbose;
+  Par.nSimulatorVerbose = nSimulatorVerbose;
+  Par.nSatSolverVerbose = nSatSolverVerbose;
+  Par.fUseBddCspf = fUseBddCspf;
+  Par.fUseBddMspf = fUseBddMspf;
+  Par.nConflictLimit = nConflictLimit;
+  Par.nSortType = nSortType;
+  Par.nOptimizerFlow = nOptimizerFlow;
+  Par.nSchedulerFlow = nSchedulerFlow;
+  rrr::Perform(&ntk, &Par);
+  Gia_Man_t *pNew = rrr::CreateGia(&ntk);
+  return pNew;
+}
+
+ABC_NAMESPACE_IMPL_END
diff --git a/src/aig/gia/module.make b/src/aig/gia/module.make
index 27240f443..2e64d1d1b 100644
--- a/src/aig/gia/module.make
+++ b/src/aig/gia/module.make
@@ -73,6 +73,7 @@ SRC +=    src/aig/gia/giaAig.c \
     src/aig/gia/giaResub6.c \
     src/aig/gia/giaRetime.c \
     src/aig/gia/giaRex.c \
+    src/aig/gia/giaRrr.cpp \
     src/aig/gia/giaSatEdge.c \
     src/aig/gia/giaSatLE.c \
     src/aig/gia/giaSatLut.c \
@@ -111,4 +112,4 @@ SRC +=    src/aig/gia/giaAig.c \
     src/aig/gia/giaTtopt.cpp \
     src/aig/gia/giaUnate.c \
     src/aig/gia/giaUtil.c \
-    src/aig/gia/giaBound.c
\ No newline at end of file
+    src/aig/gia/giaBound.c
diff --git a/src/base/abci/abc.c b/src/base/abci/abc.c
index ff1fa79fd..9d4ef58bf 100644
--- a/src/base/abci/abc.c
+++ b/src/base/abci/abc.c
@@ -521,6 +521,7 @@ static int Abc_CommandAbc9LNetOpt            ( Abc_Frame_t * pAbc, int argc, cha
 static int Abc_CommandAbc9Ttopt              ( Abc_Frame_t * pAbc, int argc, char ** argv );
 static int Abc_CommandAbc9Transduction       ( Abc_Frame_t * pAbc, int argc, char ** argv );
 static int Abc_CommandAbc9TranStoch          ( Abc_Frame_t * pAbc, int argc, char ** argv );
+static int Abc_CommandAbc9Rrr                ( Abc_Frame_t * pAbc, int argc, char ** argv );
 //#endif
 static int Abc_CommandAbc9LNetMap            ( Abc_Frame_t * pAbc, int argc, char ** argv );
 static int Abc_CommandAbc9Unmap              ( Abc_Frame_t * pAbc, int argc, char ** argv );
@@ -1323,6 +1324,7 @@ void Abc_Init( Abc_Frame_t * pAbc )
     Cmd_CommandAdd( pAbc, "ABC9",         "&ttopt",        Abc_CommandAbc9Ttopt,        0 );
     Cmd_CommandAdd( pAbc, "ABC9",         "&transduction", Abc_CommandAbc9Transduction, 0 );
     Cmd_CommandAdd( pAbc, "ABC9",         "&transtoch"   , Abc_CommandAbc9TranStoch,    0 );
+    Cmd_CommandAdd( pAbc, "ABC9",         "&rrr",          Abc_CommandAbc9Rrr,          0 );
 //#endif
     Cmd_CommandAdd( pAbc, "ABC9",         "&lnetmap",      Abc_CommandAbc9LNetMap,      0 );
     Cmd_CommandAdd( pAbc, "ABC9",         "&unmap",        Abc_CommandAbc9Unmap,        0 );
@@ -45112,6 +45114,132 @@ usage:
     return 1;
 }
 
+/**Function*************************************************************
+
+  Synopsis    []
+
+  Description []
+
+  SideEffects []
+
+  SeeAlso     []
+
+***********************************************************************/
+int Abc_CommandAbc9Rrr( Abc_Frame_t * pAbc, int argc, char ** argv )
+{
+    Gia_Man_t *pNew;
+    int c;
+    int iSeed = 0, nWords = 10, nTimeout = 0, nSchedulerVerbose = 1, nOptimizerVerbose = 0, nAnalyzerVerbose = 0, nSimulatorVerbose = 0, nSatSolverVerbose = 0, fUseBddCspf = 0, fUseBddMspf = 0, nConflictLimit = 0, nSortType = 12, nOptimizerFlow = 0, nSchedulerFlow = 0;
+    Extra_UtilGetoptReset();
+    while ( ( c = Extra_UtilGetopt( argc, argv, "XYRWTCGSOAPQabh" ) ) != EOF )
+    {
+        switch ( c )
+        {
+        case 'X':
+            nOptimizerFlow = atoi(argv[globalUtilOptind]);
+            globalUtilOptind++;
+            break;
+        case 'Y':
+            nSchedulerFlow = atoi(argv[globalUtilOptind]);
+            globalUtilOptind++;
+            break;
+        case 'R':
+            iSeed = atoi(argv[globalUtilOptind]);
+            globalUtilOptind++;
+            break;
+        case 'W':
+            nWords = atoi(argv[globalUtilOptind]);
+            globalUtilOptind++;
+            break;
+        case 'T':
+            nTimeout = atoi(argv[globalUtilOptind]);
+            globalUtilOptind++;
+            break;
+        case 'C':
+            nConflictLimit = atoi(argv[globalUtilOptind]);
+            globalUtilOptind++;
+            break;
+        case 'G':
+            nSortType = atoi(argv[globalUtilOptind]);
+            globalUtilOptind++;
+            break;
+        case 'S':
+            nSchedulerVerbose = atoi(argv[globalUtilOptind]);
+            globalUtilOptind++;
+            break;
+        case 'O':
+            nOptimizerVerbose = atoi(argv[globalUtilOptind]);
+            globalUtilOptind++;
+            break;
+        case 'A':
+            nAnalyzerVerbose = atoi(argv[globalUtilOptind]);
+            globalUtilOptind++;
+            break;
+        case 'P':
+            nSimulatorVerbose = atoi(argv[globalUtilOptind]);
+            globalUtilOptind++;
+            break;
+        case 'Q':
+            nSatSolverVerbose = atoi(argv[globalUtilOptind]);
+            globalUtilOptind++;
+            break;
+        case 'a':
+            fUseBddCspf ^= 1;
+            break;
+        case 'b':
+            fUseBddMspf ^= 1;
+            break;
+        case 'h':
+            goto usage;
+        default:
+            goto usage;
+        }
+    }
+
+    if ( argc > globalUtilOptind ) {
+        Abc_Print( -1, "Wrong number of auguments.\n" );
+        goto usage;
+    }
+
+    if ( pAbc->pGia == NULL )
+    {
+        Abc_Print( -1, "Empty GIA network.\n" );
+        return 1;
+    }
+    
+    pNew = Gia_ManRrr( pAbc->pGia, iSeed, nWords, nTimeout, nSchedulerVerbose, nOptimizerVerbose, nAnalyzerVerbose, nSimulatorVerbose, nSatSolverVerbose, fUseBddCspf, fUseBddMspf, nConflictLimit, nSortType, nOptimizerFlow, nSchedulerFlow );
+
+    Abc_FrameUpdateGia( pAbc, pNew );
+    
+    return 0;
+
+usage:
+      Abc_Print( -2, "usage: rrr [-XYRWTCGSOAPQ num] [-abh]\n" );
+      Abc_Print( -2, "\t        perform optimization\n" );
+      Abc_Print( -2, "\t-X num : method [default = %d]\n", nOptimizerFlow );
+      Abc_Print( -2, "\t                0: single-add resub\n" );
+      Abc_Print( -2, "\t                1: multi-add resub\n" );
+      Abc_Print( -2, "\t                2: repeat 0 and 1\n" );
+      Abc_Print( -2, "\t-Y num : flow [default = %d]\n", nSchedulerFlow );
+      Abc_Print( -2, "\t                0: apply method once\n" );
+      Abc_Print( -2, "\t                1: iterate like transtoch\n" );
+      Abc_Print( -2, "\t                2: iterate like deepsyn\n" );
+      Abc_Print( -2, "\t-R num : random number generator seed [default = %d]\n", iSeed );
+      Abc_Print( -2, "\t-W num : number of simulation words [default = %d]\n", nWords );
+      Abc_Print( -2, "\t-T num : timeout in seconds (0 = no timeout) [default = %d]\n", nTimeout );
+      Abc_Print( -2, "\t-C num : conflict limit [default = %d]\n", nConflictLimit );
+      Abc_Print( -2, "\t-G num : fanin cost function [default = %d]\n", nSortType );
+      Abc_Print( -2, "\t-S num : scheduler verbosity level [default = %d]\n", nSchedulerVerbose );
+      Abc_Print( -2, "\t-O num : optimizer verbosity level [default = %d]\n", nOptimizerVerbose );
+      Abc_Print( -2, "\t-A num : analyzer verbosity level [default = %d]\n", nAnalyzerVerbose );
+      Abc_Print( -2, "\t-P num : simulator verbosity level [default = %d]\n", nSimulatorVerbose );
+      Abc_Print( -2, "\t-Q num : SAT solver verbosity level [default = %d]\n", nSatSolverVerbose );
+      Abc_Print( -2, "\t-a     : use BDD-based analyzer (CSPF) [default = %s]\n", fUseBddCspf? "yes": "no" );
+      Abc_Print( -2, "\t-b     : use BDD-based analyzer (MSPF) [default = %s]\n", fUseBddMspf? "yes": "no" );
+      Abc_Print( -2, "\t-h     : print the command usage\n");
+    return 1;
+}
+
 /**Function*************************************************************
 
   Synopsis    []
diff --git a/src/opt/rrr/rrr.h b/src/opt/rrr/rrr.h
new file mode 100644
index 000000000..d306da3af
--- /dev/null
+++ b/src/opt/rrr/rrr.h
@@ -0,0 +1,37 @@
+#pragma once
+
+#include <cassert>
+
+#include "rrrParameter.h"
+#include "rrrTypes.h"
+#include "rrrAndNetwork.h"
+#include "rrrScheduler.h"
+#include "rrrOptimizer.h"
+#include "rrrAnalyzer.h"
+#include "rrrBddAnalyzer.h"
+#include "rrrBddMspfAnalyzer.h"
+#include "rrrSatSolver.h"
+#include "rrrSimulator.h"
+
+ABC_NAMESPACE_CXX_HEADER_START
+
+namespace rrr {
+
+  template <typename Ntk>
+  void Perform(Ntk *pNtk, Parameter const *pPar) {
+    assert(!pPar->fUseBddCspf || !pPar->fUseBddMspf);
+    if(pPar->fUseBddCspf) {
+      Scheduler<Ntk, rrr::Optimizer<Ntk, rrr::BddAnalyzer<Ntk>>> sch(pNtk, pPar);
+      sch.Run();
+    } else if(pPar->fUseBddMspf) {
+      Scheduler<Ntk, rrr::Optimizer<Ntk, rrr::BddMspfAnalyzer<Ntk>>> sch(pNtk, pPar);
+      sch.Run();
+    } else {
+      Scheduler<Ntk, rrr::Optimizer<Ntk, rrr::Analyzer<Ntk, rrr::Simulator<Ntk>, rrr::SatSolver<Ntk>>>> sch(pNtk, pPar);
+      sch.Run();
+    }
+  }
+  
+}
+
+ABC_NAMESPACE_CXX_HEADER_END
diff --git a/src/opt/rrr/rrrAbc.h b/src/opt/rrr/rrrAbc.h
new file mode 100644
index 000000000..5fc2e1ff5
--- /dev/null
+++ b/src/opt/rrr/rrrAbc.h
@@ -0,0 +1,80 @@
+#pragma once
+
+#include <string>
+
+#include "aig/gia/gia.h"
+#include "base/main/main.h"
+#include "base/cmd/cmd.h"
+
+ABC_NAMESPACE_CXX_HEADER_START
+
+namespace rrr {
+
+  template <typename Ntk>
+  void GiaReader(Gia_Man_t *pGia, Ntk *pNtk) {
+    int i;
+    Gia_Obj_t *pObj;
+    pNtk->Reserve(Gia_ManObjNum(pGia));
+    Gia_ManConst0(pGia)->Value = pNtk->GetConst0();
+    Gia_ManForEachObj1(pGia, pObj, i) {
+      if(Gia_ObjIsCi(pObj)) {
+        pObj->Value = pNtk->AddPi();
+      } else if(Gia_ObjIsCo(pObj)) {
+        pNtk->AddPo(Gia_ObjFanin0(pObj)->Value, Gia_ObjFaninC0(pObj));
+      } else {
+        // TODO: support XOR (and BUF and MUX?), maybe create another function
+        pObj->Value = pNtk->AddAnd(Gia_ObjFanin0(pObj)->Value, Gia_ObjFanin1(pObj)->Value, Gia_ObjFaninC0(pObj), Gia_ObjFaninC1(pObj));
+      }
+    }
+  }
+
+  template <typename Ntk>
+  Gia_Man_t *CreateGia(Ntk *pNtk) {
+    Gia_Man_t *pGia = Gia_ManStart(pNtk->GetNumNodes());
+    Gia_ManHashStart(pGia);
+    std::vector<int> v(pNtk->GetNumNodes());
+    v[0] = Gia_ManConst0Lit();
+    pNtk->ForEachPi([&](int id) {
+      v[id] = Gia_ManAppendCi(pGia);
+    });
+    pNtk->ForEachInt([&](int id) {
+      assert(pNtk->GetNodeType(id) == rrr::AND);
+      int x = -1;
+      pNtk->ForEachFanin(id, [&](int fi, bool c) {
+        if(x == -1) {
+          x = Abc_LitNotCond(v[fi], c);
+        } else {
+          x = Gia_ManHashAnd(pGia, x, Abc_LitNotCond(v[fi], c));
+        }
+      });
+      if(x == -1) {
+        x = Abc_LitNot(v[0]);
+      }
+      v[id] = x;
+    });
+    pNtk->ForEachPoDriver([&](int fi, bool c) {
+      Gia_ManAppendCo(pGia, Abc_LitNotCond(v[fi], c));
+    });
+    Gia_ManHashStop(pGia);
+    return pGia;  
+  }
+
+  template <typename Ntk>
+  void Abc9Execute(Ntk *pNtk, std::string Command) {
+    Abc_Frame_t *pAbc = Abc_FrameGetGlobalFrame();
+    Abc_FrameUpdateGia(pAbc, CreateGia(pNtk));
+    if(Abc_FrameIsBatchMode()) {
+      int r = Cmd_CommandExecute(pAbc, Command.c_str());
+      assert(r == 0);
+    } else {
+      Abc_FrameSetBatchMode(1);
+      int r = Cmd_CommandExecute(pAbc, Command.c_str());
+      assert(r == 0);
+      Abc_FrameSetBatchMode(0);
+    }
+    pNtk->Read(Abc_FrameReadGia(pAbc), GiaReader<Ntk>);
+  }
+
+}
+
+ABC_NAMESPACE_CXX_HEADER_END
diff --git a/src/opt/rrr/rrrAnalyzer.h b/src/opt/rrr/rrrAnalyzer.h
new file mode 100644
index 000000000..a01d811c7
--- /dev/null
+++ b/src/opt/rrr/rrrAnalyzer.h
@@ -0,0 +1,105 @@
+#pragma once
+
+#include <iostream>
+
+#include "rrrParameter.h"
+
+ABC_NAMESPACE_CXX_HEADER_START
+
+namespace rrr {
+
+  template <typename Ntk, typename Sim, typename Sol>
+  class Analyzer {
+  private:
+    // pointer to network
+    Ntk *pNtk;
+
+    // parameters
+    bool nVerbose;
+
+    // data
+    Sim sim;
+    Sol sol;
+
+  public:
+    // constructors
+    Analyzer(Ntk *pNtk, Parameter const *pPar);
+    void UpdateNetwork(Ntk *pNtk_, bool fSame);
+
+    // checks
+    bool CheckRedundancy(int id, int idx);
+    bool CheckFeasibility(int id, int fi, bool c);
+  };
+
+  /* {{{ Constructors */
+  
+  template <typename Ntk, typename Sim, typename Sol>
+  Analyzer<Ntk, Sim, Sol>::Analyzer(Ntk *pNtk, Parameter const *pPar) :
+    pNtk(pNtk),
+    nVerbose(pPar->nAnalyzerVerbose),
+    sim(pNtk, pPar),
+    sol(pNtk, pPar) {
+  }
+  
+  template <typename Ntk, typename Sim, typename Sol>
+  void Analyzer<Ntk, Sim, Sol>::UpdateNetwork(Ntk *pNtk_, bool fSame) {
+    pNtk = pNtk_;
+    sim.UpdateNetwork(pNtk, fSame);
+    sol.UpdateNetwork(pNtk, fSame);
+  }
+
+  /* }}} */
+
+  /* {{{ Checks */
+  
+  template <typename Ntk, typename Sim, typename Sol>
+  bool Analyzer<Ntk, Sim, Sol>::CheckRedundancy(int id, int idx) {
+    if(sim.CheckRedundancy(id, idx)) {
+      if(nVerbose) {
+        std::cout << "node " << id << " fanin " << (pNtk->GetCompl(id, idx)? "!": "") << pNtk->GetFanin(id, idx) << " index " << idx << " seems redundant" << std::endl;
+      }
+      SatResult r = sol.CheckRedundancy(id, idx);
+      if(r == UNSAT) {
+        if(nVerbose) {
+          std::cout << "node " << id << " fanin " << (pNtk->GetCompl(id, idx)? "!": "") << pNtk->GetFanin(id, idx) << " index " << idx << " is redundant" << std::endl;
+        }
+        return true;
+      }
+      if(r == SAT) {
+        if(nVerbose) {
+          std::cout << "node " << id << " fanin " << (pNtk->GetCompl(id, idx)? "!": "") << pNtk->GetFanin(id, idx) << " index " << idx << " is NOT redundant" << std::endl;
+        }
+        sim.AddCex(sol.GetCex());
+      }
+    }
+    return false;
+  }
+  
+  template <typename Ntk, typename Sim, typename Sol>
+  bool Analyzer<Ntk, Sim, Sol>::CheckFeasibility(int id, int fi, bool c) {
+    if(sim.CheckFeasibility(id, fi, c)) {
+      if(nVerbose) {
+        std::cout << "node " << id << " fanin " << (c? "!": "") << fi << " seems feasible" << std::endl;
+      }
+      SatResult r = sol.CheckFeasibility(id, fi, c);
+      if(r == UNSAT) {
+        if(nVerbose) {
+          std::cout << "node " << id << " fanin " << (c? "!": "") << fi << " is feasible" << std::endl;
+        }
+        return true;
+      }
+      if(r == SAT) {
+        if(nVerbose) {
+          std::cout << "node " << id << " fanin " << (c? "!": "") << fi << " is NOT feasible" << std::endl;
+        }
+        sim.AddCex(sol.GetCex());
+      }
+    }
+    return false;
+  }
+
+  /* }}} */
+  
+}
+
+ABC_NAMESPACE_CXX_HEADER_END
diff --git a/src/opt/rrr/rrrAndNetwork.h b/src/opt/rrr/rrrAndNetwork.h
new file mode 100644
index 000000000..2d321186c
--- /dev/null
+++ b/src/opt/rrr/rrrAndNetwork.h
@@ -0,0 +1,1132 @@
+#pragma once
+
+#include <iostream>
+#include <vector>
+#include <list>
+#include <set>
+#include <initializer_list>
+#include <string>
+#include <functional>
+#include <algorithm>
+#include <limits>
+
+#include "rrrParameter.h"
+#include "rrrTypes.h"
+
+ABC_NAMESPACE_CXX_HEADER_START
+
+namespace rrr {
+
+  class AndNetwork {
+  private:
+    // aliases
+    using itr = std::list<int>::iterator;
+    using citr = std::list<int>::const_iterator;
+    using ritr = std::list<int>::reverse_iterator;
+    using critr = std::list<int>::const_reverse_iterator;
+    using Callback = std::function<void(Action const &)>;
+
+    // network data
+    int nNodes; // number of allocated nodes
+    std::vector<int> vPis;
+    std::list<int> lInts; // internal nodes in topological order
+    std::set<int> sInts; // internal nodes as a set
+    std::vector<int> vPos;
+    std::vector<std::vector<int>> vvFaninEdges; // complementable edges, no duplicated fanins allowed (including complements), and nodes without fanins are treated as const-1
+    std::vector<int> vRefs; // reference count (number of fanouts)
+
+    // mark for network traversal
+    bool fLockTrav;
+    unsigned iTrav;
+    std::vector<unsigned> vTrav;
+
+    // flag used during constant propagation
+    bool fPropagating;
+
+    // callback functions
+    std::vector<Callback> vCallbacks;
+
+    // network backups
+    std::vector<AndNetwork> vBackups;
+
+    // conversion between node and edge
+    int  Node2Edge(int id, bool c) const { return (id << 1) + (int)c; }
+    int  Edge2Node(int f)          const { return f >> 1;             }
+    bool EdgeIsCompl(int f)        const { return f & 1;              }
+
+    // other private functions
+    int  CreateNode();
+    void SortInts(itr it);
+    void StartTraversal();
+    void EndTraversal();
+    void TakenAction(Action const &action) const;
+
+  public:
+    // constructors
+    AndNetwork();
+    AndNetwork(AndNetwork const &x);
+    AndNetwork &operator=(AndNetwork const &x);
+
+    // initialization APIs (should not be called after optimization has started)
+    void Clear();
+    void Reserve(int nReserve);
+    int  AddPi();
+    int  AddPo(int id, bool c);
+    int  AddAnd(int id0, int id1, bool c0, bool c1);
+    template <typename Ntk, typename Reader>
+    void Read(Ntk *pFrom, Reader &reader);
+
+    // network properties
+    bool UseComplementedEdges() const;
+    int  GetNumNodes() const; // number of allocated nodes (max id + 1)
+    int  GetNumPis() const;
+    int  GetNumInts() const;
+    int  GetNumPos() const;
+    int  GetConst0() const;
+    int  GetPi(int idx) const;
+    std::vector<int> GetPis() const;
+    std::vector<int> GetInts() const;
+    std::vector<int> GetPisInts() const;
+    std::vector<int> GetPos() const;
+    
+    // node properties
+    bool IsPi(int id) const;
+    bool IsInt(int id) const;
+    bool IsPo(int id) const;
+    NodeType GetNodeType(int id) const;
+    bool IsPoDriver(int id) const;
+    int  GetPiIndex(int id) const;
+    int  GetIntIndex(int id) const;
+    int  GetPoIndex(int id) const;
+    int  GetNumFanins(int id) const;
+    int  GetNumFanouts(int id) const;
+    int  GetFanin(int id, int idx) const;
+    bool GetCompl(int id, int idx) const;
+    int  FindFanin(int id, int fi) const;
+    bool IsReconvergent(int id);
+
+    // network traversal
+    void ForEachPi(std::function<void(int)> const &func) const;
+    void ForEachInt(std::function<void(int)> const &func) const;
+    void ForEachIntReverse(std::function<void(int)> const &func) const;
+    void ForEachPo(std::function<void(int)> const &func) const;
+    void ForEachPoDriver(std::function<void(int, bool)> const &func) const;
+    void ForEachFanin(int id, std::function<void(int, bool)> const &func) const;
+    void ForEachFaninIdx(int id, std::function<void(int, int, bool)> const &func) const; // func(fi, c, index of fi in fanin list of id)
+    void ForEachFanout(int id, bool fPos, std::function<void(int, bool)> const &func) const;
+    void ForEachFanoutRidx(int id, bool fPos, std::function<void(int, bool, int)> const &func) const; // func(fo, c, index of id in fanin list of fo)
+    void ForEachTfo(int id, bool fPos, std::function<void(int)> const &func);
+    void ForEachTfoReverse(int id, bool fPos, std::function<void(int)> const &func);
+    void ForEachTfoUpdate(int id, bool fPos, std::function<bool(int)> const &func);
+    template <template <typename...> typename Container, typename ... Ts>
+    void ForEachTfos(Container<Ts...> const &ids, bool fPos, std::function<void(int)> const &func);
+    template <template <typename...> typename Container, typename ... Ts>
+    void ForEachTfosUpdate(Container<Ts...> const &ids, bool fPos, std::function<bool(int)> const &func);
+
+    // Actions
+    void RemoveFanin(int id, int idx);
+    void RemoveUnused(int id, bool fRecursive = false);
+    void RemoveBuffer(int id);
+    void RemoveConst(int id);
+    void AddFanin(int id, int fi, bool c);
+    void TrivialCollapse(int id);
+    void TrivialDecompose(int id);
+    void SortFanins(int id, std::function<bool(int, bool, int, bool)> const &cost);
+
+    // Network cleanup
+    void Propagate(int id = -1); // all nodes unless specified
+    void Sweep(bool fPropagate);
+
+    // save & load
+    int  Save(int slot = -1); // slot is assigned automatically unless specified
+    void Load(int slot);
+    void PopBack(); // deletes the last entry of backups
+
+    // misc
+    void AddCallback(Callback const &callback);
+    void Print() const;
+  };
+
+  /* {{{ Private functions */
+
+  inline int AndNetwork::CreateNode() {
+    // TODO: reuse already allocated but dead nodes? or perform garbage collection?
+    vvFaninEdges.emplace_back();
+    vRefs.push_back(0);
+    assert(nNodes != std::numeric_limits<int>::max());
+    return nNodes++;
+  }
+
+  void AndNetwork::SortInts(itr it) {
+    ForEachFanin(*it, [&](int fi, bool c) {
+      itr it2 = std::find(it, lInts.end(), fi);
+      if(it2 != lInts.end()) {
+        lInts.erase(it2);
+        it2 = lInts.insert(it, fi);
+        SortInts(it2);
+      }
+    });
+  }
+
+  inline void AndNetwork::StartTraversal() {
+    assert(!fLockTrav);
+    fLockTrav = true;
+    vTrav.resize(nNodes);
+    iTrav++;
+    assert(iTrav != 0); //TODO: handle this overflow
+  }
+  
+  inline void AndNetwork::EndTraversal() {
+    assert(fLockTrav);
+    fLockTrav = false;
+  }
+  
+  inline void AndNetwork::TakenAction(Action const &action) const {
+    for(Callback const &callback: vCallbacks) {
+      callback(action);
+    }
+  }
+
+  /* }}} */
+
+  /* {{{ Constructors */
+
+  AndNetwork::AndNetwork() :
+    nNodes(0),
+    fLockTrav(false),
+    iTrav(0),
+    fPropagating(false) {
+    // add constant node
+    vvFaninEdges.emplace_back();
+    vRefs.push_back(0);
+    nNodes++;
+  }
+
+  AndNetwork::AndNetwork(AndNetwork const &x) :
+    fLockTrav(false),
+    iTrav(0),
+    fPropagating(false) {
+    nNodes       = x.nNodes;
+    vPis         = x.vPis;
+    lInts        = x.lInts;
+    sInts        = x.sInts;
+    vPos         = x.vPos;
+    vvFaninEdges = x.vvFaninEdges;
+    vRefs        = x.vRefs;
+  }
+  
+  AndNetwork &AndNetwork::operator=(AndNetwork const &x) {
+    nNodes       = x.nNodes;
+    vPis         = x.vPis;
+    lInts        = x.lInts;
+    sInts        = x.sInts;
+    vPos         = x.vPos;
+    vvFaninEdges = x.vvFaninEdges;
+    vRefs        = x.vRefs;
+    return *this;
+  }
+  
+  /* }}} */
+
+  /* {{{ Initialization APIs */
+
+  void AndNetwork::Clear() {
+    nNodes = 0;
+    vPis.clear();
+    lInts.clear();
+    sInts.clear();
+    vPos.clear();
+    vvFaninEdges.clear();
+    vRefs.clear();
+    fLockTrav = false;
+    iTrav = 0;
+    vTrav.clear();
+    fPropagating = false;
+    vCallbacks.clear();
+    vBackups.clear();
+    // add constant node
+    vvFaninEdges.emplace_back();
+    vRefs.push_back(0);
+    nNodes++;
+  }
+
+  void AndNetwork::Reserve(int nReserve) {
+    vvFaninEdges.reserve(nReserve);
+    vRefs.reserve(nReserve);
+  }
+  
+  inline int AndNetwork::AddPi() {
+    vPis.push_back(nNodes);
+    vvFaninEdges.emplace_back();
+    vRefs.push_back(0);
+    assert(nNodes != std::numeric_limits<int>::max());
+    return nNodes++;
+  }
+  
+  inline int AndNetwork::AddAnd(int id0, int id1, bool c0, bool c1) {
+    assert(id0 < nNodes);
+    assert(id1 < nNodes);
+    assert(id0 != id1);
+    lInts.push_back(nNodes);
+    sInts.insert(nNodes);
+    vRefs[id0]++;
+    vRefs[id1]++;
+    vvFaninEdges.emplace_back(std::initializer_list<int>({Node2Edge(id0, c0), Node2Edge(id1, c1)}));
+    vRefs.push_back(0);
+    assert(nNodes != std::numeric_limits<int>::max());
+    return nNodes++;
+  }
+
+  inline int AndNetwork::AddPo(int id, bool c) {
+    assert(id < nNodes);
+    vPos.push_back(nNodes);
+    vRefs[id]++;
+    vvFaninEdges.emplace_back(std::initializer_list<int>({Node2Edge(id, c)}));
+    vRefs.push_back(0);
+    assert(nNodes != std::numeric_limits<int>::max());
+    return nNodes++;
+  }
+
+  template <typename Ntk, typename Reader>
+  void AndNetwork::Read(Ntk *pFrom, Reader &reader) {
+    Clear();
+    reader(pFrom, this);
+  }
+  
+  /* }}} */
+  
+  /* {{{ Network properties */
+  
+  inline bool AndNetwork::UseComplementedEdges() const {
+    return true;
+  }
+  
+  inline int AndNetwork::GetNumNodes() const {
+    return nNodes;
+  }
+
+  inline int AndNetwork::GetNumPis() const {
+    assert(vPis.size() <= (std::vector<int>::size_type)std::numeric_limits<int>::max());
+    return vPis.size();
+  }
+  
+  inline int AndNetwork::GetNumInts() const {
+    assert(lInts.size() <= (std::vector<int>::size_type)std::numeric_limits<int>::max());
+    return lInts.size();
+  }
+  
+  inline int AndNetwork::GetNumPos() const {
+    assert(vPos.size() <= (std::vector<int>::size_type)std::numeric_limits<int>::max());
+    return vPos.size();
+  }
+
+  inline int AndNetwork::GetConst0() const {
+    return 0;
+  }
+  
+  inline int AndNetwork::GetPi(int idx) const {
+    return vPis[idx];
+  }
+
+  inline std::vector<int> AndNetwork::GetPis() const {
+    return vPis;
+  }
+
+  inline std::vector<int> AndNetwork::GetInts() const {
+    return std::vector<int>(lInts.begin(), lInts.end());
+  }
+  
+  inline std::vector<int> AndNetwork::GetPisInts() const {
+    std::vector<int> vPisInts = vPis;
+    vPisInts.insert(vPisInts.end(), lInts.begin(), lInts.end());
+    return vPisInts;
+  }
+  
+  inline std::vector<int> AndNetwork::GetPos() const {
+    return vPos;
+  }
+  
+  /* }}} */
+
+  /* {{{ Node properties */
+  
+  inline bool AndNetwork::IsPi(int id) const {
+    return GetNumFanins(id) == 0 && std::find(vPis.begin(), vPis.end(), id) != vPis.end();
+  }
+
+  inline bool AndNetwork::IsInt(int id) const {
+    return sInts.count(id);
+  }
+
+  inline bool AndNetwork::IsPo(int id) const {
+    return GetNumFanouts(id) == 0 && std::find(vPos.begin(), vPos.end(), id) != vPos.end();
+  }
+  
+  inline NodeType AndNetwork::GetNodeType(int id) const {
+    if(IsPi(id)) {
+      return PI;
+    }
+    if(IsPo(id)) {
+      return PO;
+    }
+    return AND;
+  }
+
+  inline bool AndNetwork::IsPoDriver(int id) const {
+    for(int po: vPos) {
+      if(GetFanin(po, 0) == id) {
+        return true;
+      }
+    }
+    return false;
+  }
+
+  int AndNetwork::GetPiIndex(int id) const {
+    assert(IsPi(id));
+    assert(vPis.size() <= (std::vector<int>::size_type)std::numeric_limits<int>::max());
+    std::vector<int>::const_iterator it = std::find(vPis.begin(), vPis.end(), id);
+    assert(it != vPis.end());
+    return std::distance(vPis.begin(), it);
+  }
+  
+  int AndNetwork::GetIntIndex(int id) const {
+    int index = 0;
+    citr it = lInts.begin();
+    for(; it != lInts.end(); it++) {
+      if(*it == id) {
+        break;
+      }
+      assert(index < (std::vector<int>::size_type)std::numeric_limits<int>::max());
+      index++;
+    }
+    assert(it != lInts.end());
+    return index;
+  }
+
+  int AndNetwork::GetPoIndex(int id) const {
+    assert(IsPo(id));
+    assert(vPos.size() <= (std::vector<int>::size_type)std::numeric_limits<int>::max());
+    std::vector<int>::const_iterator it = std::find(vPos.begin(), vPos.end(), id);
+    assert(it != vPos.end());
+    return std::distance(vPos.begin(), it);
+  }
+  
+  inline int AndNetwork::GetNumFanins(int id) const {
+    //assert(vvFaninEdges[id].size() <= (std::vector<int>::size_type)std::numeric_limits<int>::max());
+    return vvFaninEdges[id].size();
+  }
+
+  inline int AndNetwork::GetNumFanouts(int id) const {
+    return vRefs[id];
+  }
+
+  inline int AndNetwork::GetFanin(int id, int idx) const {
+    return Edge2Node(vvFaninEdges[id][idx]);
+  }
+
+  inline bool AndNetwork::GetCompl(int id, int idx) const {
+    return EdgeIsCompl(vvFaninEdges[id][idx]);
+  }
+
+  inline int AndNetwork::FindFanin(int id, int fi) const {
+    for(int idx = 0; idx < GetNumFanins(id); idx++) {
+      if(GetFanin(id, idx) == fi) {
+        return idx;
+      }
+    }
+    return -1;
+  }
+  
+  inline bool AndNetwork::IsReconvergent(int id) {
+    if(GetNumFanouts(id) <= 1) {
+      return false;
+    }
+    StartTraversal();
+    unsigned iTravStart = iTrav;
+    ForEachFanout(id, false, [&](int fo, bool c) {
+      vTrav[fo] = iTrav;
+      iTrav++;
+      assert(iTrav != 0); //TODO: handle this overflow
+    });
+    iTrav--;
+    if(iTrav <= iTravStart) {
+      // less than one fanouts excluding POs
+      EndTraversal();
+      return false;
+    }
+    citr it = lInts.begin();
+    while(vTrav[*it] < iTravStart && it != lInts.end()) {
+      it++;
+    }
+    it++;
+    for(; it != lInts.end(); it++) {
+      for(int fi_edge: vvFaninEdges[*it]) {
+        int fi = Edge2Node(fi_edge);
+        if(vTrav[fi] >= iTravStart) {
+          if(vTrav[*it] >= iTravStart && vTrav[*it] != vTrav[fi]) {
+            EndTraversal();
+            return true;
+          }
+          vTrav[*it] = vTrav[fi];
+        }
+      }
+    }
+    EndTraversal();
+    return false;
+  }
+
+  /* }}} */
+
+  /* {{{ Network traversal */
+
+  inline void AndNetwork::ForEachPi(std::function<void(int)> const &func) const {
+    for(int pi: vPis) {
+      func(pi);
+    }
+  }
+  
+  inline void AndNetwork::ForEachInt(std::function<void(int)> const &func) const {
+    for(int id: lInts) {
+      func(id);
+    }
+  }
+
+  inline void AndNetwork::ForEachIntReverse(std::function<void(int)> const &func) const {
+    for(critr it = lInts.rbegin(); it != lInts.rend(); it++) {
+      func(*it);
+    }
+  }
+
+  inline void AndNetwork::ForEachPo(std::function<void(int)> const &func) const {
+    for(int po: vPos) {
+      func(po);
+    }
+  }
+  
+  inline void AndNetwork::ForEachPoDriver(std::function<void(int, bool)> const &func) const {
+    for(int po: vPos) {
+      func(GetFanin(po, 0), GetCompl(po, 0));
+    }
+  }
+  
+  inline void AndNetwork::ForEachFanin(int id, std::function<void(int, bool)> const &func) const {
+    for(int fi_edge: vvFaninEdges[id]) {
+      func(Edge2Node(fi_edge), EdgeIsCompl(fi_edge));
+    }
+  }
+
+  inline void AndNetwork::ForEachFaninIdx(int id, std::function<void(int, int, bool)> const &func) const {
+    for(int idx = 0; idx < GetNumFanins(id); idx++) {
+      func(idx, GetFanin(id, idx), GetCompl(id, idx));
+    }
+  }
+  
+  inline void AndNetwork::ForEachFanout(int id, bool fPos, std::function<void(int, bool)> const &func) const {
+    if(vRefs[id] == 0) {
+      return;
+    }
+    citr it = std::find(lInts.begin(), lInts.end(), id);
+    assert(it != lInts.end());
+    it++;
+    int nRefs = vRefs[id];
+    for(; nRefs != 0 && it != lInts.end(); it++) {
+      int idx = FindFanin(*it, id);
+      if(idx >= 0) {
+        func(*it, GetCompl(*it, idx));
+        nRefs--;
+      }
+    }
+    if(fPos && nRefs != 0) {
+      for(int po: vPos) {
+        if(GetFanin(po, 0) == id) {
+          func(po, GetCompl(po, 0));
+          nRefs--;
+          if(nRefs == 0) {
+            break;
+          }
+        }
+      }
+    }
+    assert(!fPos || nRefs == 0);
+  }
+  
+  inline void AndNetwork::ForEachFanoutRidx(int id, bool fPos, std::function<void(int, bool, int)> const &func) const {
+    if(vRefs[id] == 0) {
+      return;
+    }
+    citr it = std::find(lInts.begin(), lInts.end(), id);
+    assert(it != lInts.end());
+    it++;
+    int nRefs = vRefs[id];
+    for(; nRefs != 0 && it != lInts.end(); it++) {
+      int idx = FindFanin(*it, id);
+      if(idx >= 0) {
+        func(*it, GetCompl(*it, idx), idx);
+        nRefs--;
+      }
+    }
+    if(fPos && nRefs != 0) {
+      for(int po: vPos) {
+        if(GetFanin(po, 0) == id) {
+          func(po, GetCompl(po, 0), 0);
+          nRefs--;
+          if(nRefs == 0) {
+            break;
+          }
+        }
+      }
+    }
+    assert(!fPos || nRefs == 0);
+  }
+  
+  inline void AndNetwork::ForEachTfo(int id, bool fPos, std::function<void(int)> const &func) {
+    // this does not include id itself
+    if(vRefs[id] == 0) {
+      return;
+    }
+    StartTraversal();
+    vTrav[id] = iTrav;
+    citr it = std::find(lInts.begin(), lInts.end(), id);
+    assert(it != lInts.end());
+    it++;
+    for(; it != lInts.end(); it++) {
+      for(int fi_edge: vvFaninEdges[*it]) {
+        if(vTrav[Edge2Node(fi_edge)] == iTrav) {
+          func(*it);
+          vTrav[*it] = iTrav;
+          break;
+        }
+      }
+    }
+    if(fPos) {
+      for(int po: vPos) {
+        if(vTrav[GetFanin(po, 0)] == iTrav) {
+          func(po);
+          vTrav[po] = iTrav;
+        }
+      }
+    }
+    EndTraversal();
+  }
+
+  inline void AndNetwork::ForEachTfoReverse(int id, bool fPos, std::function<void(int)> const &func) {
+    // this does not include id itself
+    if(vRefs[id] == 0) {
+      return;
+    }
+    StartTraversal();
+    vTrav[id] = iTrav;
+    citr it = std::find(lInts.begin(), lInts.end(), id);
+    assert(it != lInts.end());
+    it++;
+    for(; it != lInts.end(); it++) {
+      for(int fi_edge: vvFaninEdges[*it]) {
+        if(vTrav[Edge2Node(fi_edge)] == iTrav) {
+          vTrav[*it] = iTrav;
+          break;
+        }
+      }
+    }
+    if(fPos) {
+      for(int po: vPos) {
+        if(vTrav[GetFanin(po, 0)] == iTrav) {
+          vTrav[po] = iTrav;
+        }
+      }
+    }
+    EndTraversal(); // release here so func can call IsReconvergent
+    unsigned iTravTfo = iTrav;
+    if(fPos) {
+      // use reverse order even for POs
+      for(std::vector<int>::const_reverse_iterator it = vPos.rbegin(); it != vPos.rend(); it++) {
+        assert(vTrav[*it] <= iTravTfo); // make sure func does not touch vTrav of preceding nodes
+        if(vTrav[*it] == iTravTfo) {
+          func(*it);
+        }
+      }
+    }
+    for(critr it = lInts.rbegin(); *it != id; it++) {
+      assert(vTrav[*it] <= iTravTfo); // make sure func does not touch vTrav of preceding nodes
+      if(vTrav[*it] == iTravTfo) {
+        func(*it);
+      }
+    }
+  }
+
+  inline void AndNetwork::ForEachTfoUpdate(int id, bool fPos, std::function<bool(int)> const &func) {
+    // this does not include id itself
+    if(vRefs[id] == 0) {
+      return;
+    }
+    StartTraversal();
+    vTrav[id] = iTrav;
+    citr it = std::find(lInts.begin(), lInts.end(), id);
+    assert(it != lInts.end());
+    it++;
+    for(; it != lInts.end(); it++) {
+      for(int fi_edge: vvFaninEdges[*it]) {
+        if(vTrav[Edge2Node(fi_edge)] == iTrav) {
+          if(func(*it)) {
+            vTrav[*it] = iTrav;
+          }
+          break;
+        }
+      }
+    }
+    if(fPos) {
+      for(int po: vPos) {
+        if(vTrav[GetFanin(po, 0)] == iTrav) {
+          if(func(po)) {
+            vTrav[po] = iTrav;
+          }
+        }
+      }
+    }
+    EndTraversal();
+  }
+
+  template <template <typename...> typename Container, typename ... Ts>
+  inline void AndNetwork::ForEachTfos(Container<Ts...> const &ids, bool fPos, std::function<void(int)> const &func) {
+    // this includes ids themselves
+    StartTraversal();
+    for(int id: ids) {
+      vTrav[id] = iTrav;
+    }
+    citr it = lInts.begin();
+    while(vTrav[*it] != iTrav && it != lInts.end()) {
+      it++;
+    }
+    for(; it != lInts.end(); it++) {
+      if(vTrav[*it] == iTrav) {
+        func(*it);
+      } else {
+        for(int fi_edge: vvFaninEdges[*it]) {
+          if(vTrav[Edge2Node(fi_edge)] == iTrav) {
+            func(*it);
+            vTrav[*it] = iTrav;
+            break;
+          }
+        }
+      }
+    }
+    if(fPos) {
+      for(int po: vPos) {
+        if(vTrav[po] == iTrav || vTrav[GetFanin(po, 0)] == iTrav) {
+          func(po);
+          vTrav[po] = iTrav;
+        }
+      }
+    }
+    EndTraversal();
+  }
+  
+  template <template <typename...> typename Container, typename ... Ts>
+  inline void AndNetwork::ForEachTfosUpdate(Container<Ts...> const &ids, bool fPos, std::function<bool(int)> const &func) {
+    // this includes ids themselves
+    StartTraversal();
+    for(int id: ids) {
+      vTrav[id] = iTrav;
+    }
+    citr it = lInts.begin();
+    while(vTrav[*it] != iTrav && it != lInts.end()) {
+      it++;
+    }
+    for(; it != lInts.end(); it++) {
+      if(vTrav[*it] == iTrav) {
+        if(!func(*it)) {
+          vTrav[*it] = 0;
+        }
+      } else {
+        for(int fi_edge: vvFaninEdges[*it]) {
+          if(vTrav[Edge2Node(fi_edge)] == iTrav) {
+            if(func(*it)) {
+              vTrav[*it] = iTrav;
+            }
+            break;
+          }
+        }
+      }
+    }
+    if(fPos) {
+      for(int po: vPos) {
+        if(vTrav[po] == iTrav) {
+          if(!func(po)) {
+            vTrav[po] = 0;
+          }
+        } else if(vTrav[GetFanin(po, 0)] == iTrav) {
+          if(func(po)) {
+            vTrav[po] = iTrav;
+          }
+        }
+      }
+    }
+    EndTraversal();
+  }
+
+  /* }}} */
+  
+  /* {{{ Actions */
+  
+  void AndNetwork::RemoveFanin(int id, int idx) {
+    Action action;
+    action.type = REMOVE_FANIN;
+    action.id = id;
+    action.idx = idx;
+    int fi = GetFanin(id, idx);
+    bool c = GetCompl(id, idx);
+    action.fi = fi;
+    action.c = c;
+    vRefs[fi]--;
+    vvFaninEdges[id].erase(vvFaninEdges[id].begin() + idx);
+    TakenAction(action);
+  }
+
+  void AndNetwork::RemoveUnused(int id, bool fRecursive) {
+    assert(vRefs[id] == 0);
+    Action action;
+    action.type = REMOVE_UNUSED;
+    action.id = id;
+    ForEachFanin(id, [&](int fi, bool c) {
+      action.vFanins.push_back(fi);
+      vRefs[fi]--;
+    });
+    vvFaninEdges[id].clear();
+    itr it = std::find(lInts.begin(), lInts.end(), id);
+    lInts.erase(it);
+    sInts.erase(id);
+    TakenAction(action);
+    if(fRecursive) {
+      for(int fi: action.vFanins) {
+        if(vRefs[fi] == 0) {
+          RemoveUnused(fi, true);
+        }
+      }
+    }
+  }
+
+  void AndNetwork::RemoveBuffer(int id) {
+    assert(GetNumFanins(id) == 1);
+    assert(!fPropagating || fLockTrav);
+    int fi = GetFanin(id, 0);
+    bool c = GetCompl(id, 0);
+    // check if it is buffering constant
+    if(fi == GetConst0()) {
+      RemoveConst(id);
+      return;
+    }
+    // remove if substitution would lead to duplication with the same polarity
+    ForEachFanoutRidx(id, false, [&](int fo, bool foc, int idx) {
+      int idx2 = FindFanin(fo, fi);
+      if(idx2 != -1 && GetCompl(fo, idx2) == (c ^ foc)) {
+        RemoveFanin(fo, idx);
+        if(fPropagating && GetNumFanins(fo) == 1) {
+          vTrav[fo] = iTrav;
+        }
+      }
+    });
+    // substitute node with fanin or const-0
+    Action action;
+    action.type = REMOVE_BUFFER;
+    action.id = id;
+    action.fi = fi;
+    action.c = c;
+    ForEachFanoutRidx(id, true, [&](int fo, bool foc, int idx) {
+      action.vFanouts.push_back(fo);
+      int idx2 = FindFanin(fo, fi);
+      if(idx2 != -1) { // substitute with const-0 in case of duplication
+        assert(GetCompl(fo, idx2) != (c ^ foc)); // of a different polarity
+        vRefs[GetConst0()]++;
+        vvFaninEdges[fo][idx] = Node2Edge(GetConst0(), 0);
+        if(fPropagating) {
+          vTrav[fo] = iTrav;
+        }
+      } else { // otherwise, substitute with fanin
+        vvFaninEdges[fo][idx] = Node2Edge(fi, c ^ foc);
+        vRefs[fi]++;
+      }
+    });
+    // remove node
+    vRefs[id] = 0;
+    vRefs[fi]--;
+    vvFaninEdges[id].clear();
+    if(!fPropagating) {
+      itr it = std::find(lInts.begin(), lInts.end(), id);
+      lInts.erase(it);
+    }
+    sInts.erase(id);
+    TakenAction(action);
+  }
+
+  void AndNetwork::RemoveConst(int id) {
+    assert(GetNumFanins(id) == 0 || FindFanin(id, GetConst0()) != -1);
+    assert(!fPropagating || fLockTrav);
+    bool c = (GetNumFanins(id) == 0);
+    // just remove immediately if polarity is true but not PO
+    ForEachFanoutRidx(id, false, [&](int fo, bool foc, int idx) {
+      if(c ^ foc) {
+        assert(!IsPo(fo));
+        RemoveFanin(fo, idx);
+        if(fPropagating && GetNumFanins(fo) <= 1) {
+          vTrav[fo] = iTrav;
+        }
+      }
+    });
+    // substitute with constant
+    Action action;
+    action.type = REMOVE_CONST;
+    action.id = id;
+    ForEachFanoutRidx(id, true, [&](int fo, bool foc, int idx) {
+      action.vFanouts.push_back(fo);
+      vRefs[GetConst0()]++;
+      vvFaninEdges[fo][idx] = Node2Edge(GetConst0(), c ^ foc);
+      if(fPropagating) {
+        vTrav[fo] = iTrav;
+      }
+    });
+    // remove node
+    vRefs[id] = 0;
+    ForEachFanin(id, [&](int fi, bool c) {
+      vRefs[fi]--;
+      action.vFanins.push_back(fi);
+    });
+    vvFaninEdges[id].clear();
+    if(!fPropagating) {
+      itr it = std::find(lInts.begin(), lInts.end(), id);
+      lInts.erase(it);
+    }
+    sInts.erase(id);
+    TakenAction(action);
+  }
+
+  void AndNetwork::AddFanin(int id, int fi, bool c) {
+    assert(FindFanin(id, fi) == -1); // no duplication
+    assert(fi != GetConst0() || !c); // no const-1
+    Action action;
+    action.type = ADD_FANIN;
+    action.id = id;
+    action.idx = vvFaninEdges[id].size();
+    action.fi = fi;
+    action.c = c;
+    itr it = std::find(lInts.begin(), lInts.end(), id);
+    itr it2 = std::find(it, lInts.end(), fi);
+    if(it2 != lInts.end()) {
+      lInts.erase(it2);
+      it2 = lInts.insert(it, fi);
+      SortInts(it2);
+    }
+    vRefs[fi]++;
+    vvFaninEdges[id].push_back(Node2Edge(fi, c));
+    TakenAction(action);
+  }
+
+  void AndNetwork::TrivialCollapse(int id) {
+    for(int idx = 0; idx < GetNumFanins(id);) {
+      int fi_edge = vvFaninEdges[id][idx];
+      int fi = Edge2Node(fi_edge);
+      bool c = EdgeIsCompl(fi_edge);
+      if(!IsPi(fi) && !c && vRefs[fi] == 1) {
+        Action action;
+        action.type = TRIVIAL_COLLAPSE;
+        action.id = id;
+        action.idx = idx;
+        action.fi = fi;
+        action.c = c;
+        std::vector<int>::iterator it = vvFaninEdges[id].begin() + idx;
+        it = vvFaninEdges[id].erase(it);
+        vvFaninEdges[id].insert(it, vvFaninEdges[fi].begin(), vvFaninEdges[fi].end());
+        ForEachFanin(fi, [&](int fi, bool c) {
+          action.vFanins.push_back(fi);
+        });
+        // remove collapsed fanin
+        vRefs[fi] = 0;
+        vvFaninEdges[fi].clear();
+        lInts.erase(std::find(lInts.begin(), lInts.end(), fi));
+        sInts.erase(fi);
+        TakenAction(action);
+      } else {
+        idx++;
+      }
+    }
+  }
+
+  void AndNetwork::TrivialDecompose(int id) {
+    while(GetNumFanins(id) > 2) {
+      Action action;
+      action.type = TRIVIAL_DECOMPOSE;
+      action.id = id;
+      action.idx = vvFaninEdges[id].size() - 2;
+      int new_fi = CreateNode();
+      action.fi = new_fi;
+      int fi_edge1 = vvFaninEdges[id].back();
+      vvFaninEdges[id].pop_back();
+      int fi_edge0 = vvFaninEdges[id].back();
+      vvFaninEdges[id].pop_back();
+      vvFaninEdges[new_fi].push_back(fi_edge0);
+      action.vFanins.push_back(Edge2Node(fi_edge0));
+      vvFaninEdges[new_fi].push_back(fi_edge1);
+      action.vFanins.push_back(Edge2Node(fi_edge1));
+      vvFaninEdges[id].push_back(Node2Edge(new_fi, false));
+      vRefs[new_fi]++;
+      itr it = std::find(lInts.begin(), lInts.end(), id);
+      lInts.insert(it, new_fi);
+      sInts.insert(new_fi);
+      TakenAction(action);
+    }
+  }
+
+  void AndNetwork::SortFanins(int id, std::function<bool(int, bool, int, bool)> const &comp) {
+    std::vector<int> vFaninEdges = vvFaninEdges[id];
+    std::sort(vvFaninEdges[id].begin(), vvFaninEdges[id].end(), [&](int i, int j) {
+      return comp(Edge2Node(i), EdgeIsCompl(i), Edge2Node(j), EdgeIsCompl(j));
+    });
+    if(vFaninEdges == vvFaninEdges[id]) {
+      return;
+    }
+    Action action;
+    action.type = SORT_FANINS;
+    action.id = id;
+    assert(vFaninEdges.size() <= (std::vector<int>::size_type)std::numeric_limits<int>::max());
+    for(int fanin_edge: vvFaninEdges[id]) {
+      std::vector<int>::const_iterator it = std::find(vFaninEdges.begin(), vFaninEdges.end(), fanin_edge);
+      action.vIndices.push_back(std::distance(vFaninEdges.cbegin(), it));
+    }
+    TakenAction(action);
+  }
+
+  /* }}} */
+
+  /* {{{ Network cleanup */
+  
+  void AndNetwork::Propagate(int id) {
+    StartTraversal();
+    itr it;
+    if(id == -1) {
+      ForEachInt([&](int id) {
+        if(GetNumFanins(id) <= 1 || FindFanin(id, GetConst0()) != -1) {
+          vTrav[id] = iTrav;
+        }
+      });
+      it = lInts.begin();
+      while(vTrav[*it] != iTrav && it != lInts.end()) {
+        it++;
+      }
+    } else {
+      vTrav[id] = iTrav;
+      it = std::find(lInts.begin(), lInts.end(), id);
+    }
+    fPropagating = true;
+    while(it != lInts.end()) {
+      if(vTrav[*it] == iTrav) {
+        if(GetNumFanins(*it) == 1) {
+          RemoveBuffer(*it);
+        } else {
+          RemoveConst(*it);
+        }
+        it = lInts.erase(it);
+      } else {
+        it++;
+      }
+    }
+    fPropagating = false;
+    EndTraversal();
+  }
+
+  void AndNetwork::Sweep(bool fPropagate) {
+    if(fPropagate) {
+      Propagate();
+    }
+    for(ritr it = lInts.rbegin(); it != lInts.rend();) {
+      if(vRefs[*it] == 0) {
+        RemoveUnused(*it);
+        it = ritr(lInts.erase(--it.base()));
+      } else {
+        it++;
+      }
+    }
+  }
+  
+  /* }}} */
+
+  /* {{{ Save & load */
+
+  int AndNetwork::Save(int slot) {
+    Action action;
+    action.type = SAVE;
+    if(slot < 0) {
+      slot = vBackups.size();
+      vBackups.push_back(*this);
+    } else {
+      assert(slot < vBackups.size());
+      vBackups[slot] = *this;
+    }
+    action.idx = slot;
+    TakenAction(action);
+    return slot;
+  }
+
+  void AndNetwork::Load(int slot) {
+    assert(slot < vBackups.size());
+    Action action;
+    action.type = LOAD;
+    action.idx = slot;
+    *this = vBackups[slot];
+    TakenAction(action);
+  }
+
+  void AndNetwork::PopBack() {
+    assert(!vBackups.empty());
+    Action action;
+    action.type = POP_BACK;
+    action.idx = vBackups.size() - 1;
+    vBackups.pop_back();
+    TakenAction(action);
+  }
+
+  /* }}} */
+  
+  /* {{{ Misc */
+
+  void AndNetwork::AddCallback(Callback const &callback) {
+    vCallbacks.push_back(callback);
+  }
+  
+  void AndNetwork::Print() const {
+    std::cout << "pi: ";
+    std::string delim = "";
+    ForEachPi([&](int id) {
+      std::cout << delim << id;
+      delim = ", ";
+    });
+    std::cout << std::endl;
+    ForEachInt([&](int id) {
+      std::cout << "node " << id << ": ";
+      delim = "";
+      ForEachFanin(id, [&](int fi, bool c) {
+        std::cout << delim;
+        if(c) {
+          std::cout << "!";
+        }
+        std::cout << fi;
+        delim = ", ";
+      });
+      std::cout << " (ref = " << vRefs[id] << ")";
+      std::cout << std::endl;
+    });
+    std::cout << "po: ";
+    delim = "";
+    ForEachPoDriver([&](int fi, bool c) {
+      std::cout << delim;
+      if(c) {
+        std::cout << "!";
+      }
+      std::cout << fi;
+      delim = ", ";
+    });
+    std::cout << std::endl;
+  }
+
+  /* }}} Misc end */
+
+}
+
+ABC_NAMESPACE_CXX_HEADER_END
diff --git a/src/opt/rrr/rrrBddAnalyzer.h b/src/opt/rrr/rrrBddAnalyzer.h
new file mode 100644
index 000000000..842e7714e
--- /dev/null
+++ b/src/opt/rrr/rrrBddAnalyzer.h
@@ -0,0 +1,667 @@
+#pragma once
+
+#include <iostream>
+#include <vector>
+
+#include <aig/gia/giaNewBdd.h>
+
+#include "rrrParameter.h"
+#include "rrrTypes.h"
+
+ABC_NAMESPACE_CXX_HEADER_START
+
+namespace rrr {
+
+  template <typename Ntk>
+  class BddAnalyzer {
+  private:
+    // aliases
+    using lit = int;
+    static constexpr lit LitMax = 0xffffffff;
+
+    // pointer to network
+    Ntk *pNtk;
+    
+    // parameters
+    int nVerbose;
+    
+    // data
+    NewBdd::Man *pBdd;
+    int target;
+    std::vector<lit> vFs;
+    std::vector<lit> vGs;
+    std::vector<std::vector<lit>> vvCs;
+    bool fResim;
+    std::vector<bool> vUpdates;
+    std::vector<bool> vGUpdates;
+    std::vector<bool> vCUpdates;
+
+    // backups
+    std::vector<BddAnalyzer> vBackups;
+
+    // BDD utils
+    void IncRef(lit x) const;
+    void DecRef(lit x) const;
+    void Assign(lit &x, lit y) const;
+    void CopyVec(std::vector<lit> &x, std::vector<lit> const &y) const;
+    void DelVec(std::vector<lit> &v) const;
+    void CopyVecVec(std::vector<std::vector<lit>> &x, std::vector<std::vector<lit>> const &y) const;
+    void DelVecVec(std::vector<std::vector<lit>> &v) const;
+    lit  Xor(lit x, lit y) const;
+
+    // callback
+    void ActionCallback(Action const &action);
+
+    // allocation
+    void Allocate();
+
+    // simulation
+    void SimulateNode(int id, std::vector<lit> &v) const;
+    void Simulate();
+    
+    // CSPF computation
+    bool ComputeG(int id);
+    void ComputeC(int id);
+    void CspfNode(int id);
+    void Cspf(int id = -1, bool fC = true);
+
+    // save & load
+    void Save(int slot);
+    void Load(int slot);
+    void PopBack();
+
+  public:
+    // constructors
+    BddAnalyzer();
+    BddAnalyzer(Ntk *pNtk, Parameter const *pPar);
+    ~BddAnalyzer();
+    void UpdateNetwork(Ntk *pNtk_, bool fSame);
+
+    // checks
+    bool CheckRedundancy(int id, int idx);
+    bool CheckFeasibility(int id, int fi, bool c);
+  };
+  
+  /* {{{ BDD utils */
+  
+  template <typename Ntk>
+  inline void BddAnalyzer<Ntk>::IncRef(lit x) const {
+    if(x != LitMax) {
+      pBdd->IncRef(x);
+    }
+  }
+  
+  template <typename Ntk>
+  inline void BddAnalyzer<Ntk>::DecRef(lit x) const {
+    if(x != LitMax) {
+      pBdd->DecRef(x);
+    }
+  }
+  
+  template <typename Ntk>
+  inline void BddAnalyzer<Ntk>::Assign(lit &x, lit y) const {
+    DecRef(x);
+    x = y;
+    IncRef(x);
+  }
+
+  template <typename Ntk>
+  inline void BddAnalyzer<Ntk>::CopyVec(std::vector<lit> &x, std::vector<lit> const &y) const {
+    for(size_t i = y.size(); i < x.size(); i++) {
+      DecRef(x[i]);
+    }
+    x.resize(y.size(), LitMax);
+    for(size_t i = 0; i < y.size(); i++) {
+      if(x[i] != y[i]) {
+        DecRef(x[i]);
+        x[i] = y[i];
+        IncRef(x[i]);
+      }
+    }
+  }
+
+  template <typename Ntk>
+  inline void BddAnalyzer<Ntk>::DelVec(std::vector<lit> &v) const {
+    for(lit x: v) {
+      DecRef(x);
+    }
+    v.clear();
+  }
+
+  template <typename Ntk>
+  inline void BddAnalyzer<Ntk>::CopyVecVec(std::vector<std::vector<lit>> &x, std::vector<std::vector<lit>> const &y) const {
+    for(size_t i = y.size(); i < x.size(); i++) {
+      DelVec(x[i]);
+    }
+    x.resize(y.size());
+    for(size_t i = 0; i < y.size(); i++) {
+      CopyVec(x[i], y[i]);
+    }
+  }
+
+  template <typename Ntk>
+  inline void BddAnalyzer<Ntk>::DelVecVec(std::vector<std::vector<lit>> &v) const {
+    for(size_t i = 0; i < v.size(); i++) {
+      DelVec(v[i]);
+    }
+    v.clear();
+  }
+  
+  template <typename Ntk>
+  inline int BddAnalyzer<Ntk>::Xor(lit x, lit y) const {
+    lit f = pBdd->And(x, pBdd->LitNot(y));
+    pBdd->IncRef(f);
+    lit g = pBdd->And(pBdd->LitNot(x), y);
+    pBdd->IncRef(g);
+    lit r = pBdd->Or(f, g);
+    pBdd->DecRef(f);
+    pBdd->DecRef(g);
+    return r;
+  }
+  
+  /* }}} */
+
+  /* {{{ Callback */
+
+  template <typename Ntk>
+  void BddAnalyzer<Ntk>::ActionCallback(Action const &action) {
+    switch(action.type) {
+    case REMOVE_FANIN:
+      vUpdates[action.id] = true;
+      vGUpdates[action.fi] = true;
+      DecRef(vvCs[action.id][action.idx]);
+      vvCs[action.id].erase(vvCs[action.id].begin() + action.idx);
+      break;
+    case REMOVE_UNUSED:
+      if(vGUpdates[action.id] || vCUpdates[action.id]) {
+        for(int fi: action.vFanins) {
+          vGUpdates[fi] = true;
+        }
+      }
+      Assign(vFs[action.id], LitMax);
+      Assign(vGs[action.id], LitMax);
+      DelVec(vvCs[action.id]);
+      break;
+    case REMOVE_BUFFER:
+      if(vUpdates[action.id]) {
+        for(int fo: action.vFanouts) {
+          vUpdates[fo] = true;
+          vCUpdates[fo] = true;
+        }
+      }
+      if(vGUpdates[action.id] || vCUpdates[action.id]) {
+        vGUpdates[action.fi] = true;
+      }
+      Assign(vFs[action.id], LitMax);
+      Assign(vGs[action.id], LitMax);
+      DelVec(vvCs[action.id]);
+      break;
+    case REMOVE_CONST:
+      if(vUpdates[action.id]) {
+        for(int fo: action.vFanouts) {
+          vUpdates[fo] = true;
+          vCUpdates[fo] = true;
+        }
+      }
+      for(int fi: action.vFanins) {
+        vGUpdates[fi] = true;
+      }
+      Assign(vFs[action.id], LitMax);
+      Assign(vGs[action.id], LitMax);
+      DelVec(vvCs[action.id]);
+      break;
+    case ADD_FANIN:
+      vUpdates[action.id] = true;
+      vCUpdates[action.id] = true;
+      vvCs[action.id].insert(vvCs[action.id].begin() + action.idx, LitMax);
+      break;
+    case TRIVIAL_COLLAPSE: {
+      if(vGUpdates[action.fi] || vCUpdates[action.fi]) {
+        vCUpdates[action.id] = true;
+      }
+      std::vector<lit>::iterator it = vvCs[action.id].begin() + action.idx;
+      DecRef(*it);
+      it = vvCs[action.id].erase(it);
+      vvCs[action.id].insert(it,  vvCs[action.fi].begin(), vvCs[action.fi].end());
+      vvCs[action.fi].clear();
+      Assign(vFs[action.fi], LitMax);
+      Assign(vGs[action.fi], LitMax);
+      break;
+    }
+    case TRIVIAL_DECOMPOSE: {
+      Allocate();
+      SimulateNode(action.fi, vFs);
+      assert(vGs[action.fi] == LitMax);
+      Assign(vGs[action.fi], vGs[action.id]);
+      std::vector<lit>::iterator it = vvCs[action.id].begin() + action.idx;
+      assert(vvCs[action.fi].empty());
+      vvCs[action.fi].insert(vvCs[action.fi].begin(), it, vvCs[action.id].end());
+      vvCs[action.id].erase(it, vvCs[action.id].end());
+      assert(vvCs[action.id].size() == action.idx);
+      vvCs[action.id].resize(action.idx + 1, LitMax);
+      Assign(vvCs[action.id][action.idx], vGs[action.fi]);
+      vUpdates[action.fi] = false;
+      vGUpdates[action.fi] = false;
+      vCUpdates[action.fi] = vCUpdates[action.id];
+      break;
+    }
+    case SORT_FANINS: {
+      std::vector<lit> vCs = vvCs[action.id];
+      vvCs[action.id].clear();
+      for(int index: action.vIndices) {
+        vvCs[action.id].push_back(vCs[index]);
+      }
+      if(!fResim && target != -1 && target != action.id) {
+        pNtk->ForEachTfo(target, false, [&](int fo) {
+          if(fResim) {
+            return;
+          }
+          if(fo == action.id) {
+            fResim = true;
+          }
+        });
+      }
+      vCUpdates[action.id] = true;
+      break;
+    }
+    case SAVE:
+      Save(action.idx);
+      break;
+    case LOAD:
+      Load(action.idx);
+      break;
+    case POP_BACK:
+      PopBack();
+      break;
+    default:
+      assert(0);
+    }
+  }
+  
+  /* }}} */
+
+  /* {{{ Allocation */
+
+  template <typename Ntk>
+  void BddAnalyzer<Ntk>::Allocate() {
+    int nNodes = pNtk->GetNumNodes();
+    vFs.resize(nNodes, LitMax);
+    vGs.resize(nNodes, LitMax);
+    vvCs.resize(nNodes);
+    vUpdates.resize(nNodes);
+    vGUpdates.resize(nNodes);
+    vCUpdates.resize(nNodes);
+  }
+
+  /* }}} */
+  
+  /* {{{ Simulation */
+  
+  template <typename Ntk>
+  inline void BddAnalyzer<Ntk>::SimulateNode(int id, std::vector<lit> &v) const {
+    if(nVerbose) {
+      std::cout << "simulating node " << id << std::endl;
+    }
+    Assign(v[id], pBdd->Const1());
+    pNtk->ForEachFanin(id, [&](int fi, bool c) {
+      Assign(v[id], pBdd->And(v[id], pBdd->LitNotCond(v[fi], c)));
+    });
+  }
+  
+  template <typename Ntk>
+  void BddAnalyzer<Ntk>::Simulate() {
+    if(nVerbose) {
+      std::cout << "symbolic simulation with BDD" << std::endl;
+    }
+    pNtk->ForEachInt([&](int id) {
+      if(vUpdates[id]) {
+        lit x = vFs[id];
+        IncRef(x);
+        SimulateNode(id, vFs);
+        DecRef(x);
+        if(!pBdd->LitIsEq(x, vFs[id])) {
+          pNtk->ForEachFanout(id, false, [&](int fo, bool c) {
+            vUpdates[fo] = true;
+            vCUpdates[fo] = true;
+          });
+        }
+        vUpdates[id] = false;
+      }
+    });
+  }
+
+  /* }}} */
+
+  /* {{{ CSPF computation */
+  
+  template <typename Ntk>
+  inline bool BddAnalyzer<Ntk>::ComputeG(int id) {
+    if(pNtk->GetNumFanouts(id) == 0) {
+      if(pBdd->IsConst1(vGs[id])) {
+        return false;
+      }
+      Assign(vGs[id], pBdd->Const1());
+      return true;
+    }
+    lit x = pBdd->Const1();
+    IncRef(x);
+    pNtk->ForEachFanoutRidx(id, true, [&](int fo, bool c, int idx) {
+      Assign(x, pBdd->And(x, vvCs[fo][idx]));
+    });
+    if(pBdd->LitIsEq(vGs[id], x)) {
+      DecRef(x);
+      return false;
+    }
+    Assign(vGs[id], x);
+    DecRef(x);
+    return true;
+  }
+
+  template <typename Ntk>
+  inline void BddAnalyzer<Ntk>::ComputeC(int id) {
+    int nFanins = pNtk->GetNumFanins(id);
+    assert(vvCs[id].size() == nFanins);
+    if(pBdd->IsConst1(vGs[id])) {
+      for(int idx = 0; idx < nFanins; idx++) {
+        if(!pBdd->IsConst1(vvCs[id][idx])) {
+          Assign(vvCs[id][idx], pBdd->Const1());
+          int fi = pNtk->GetFanin(id, idx);
+          vGUpdates[fi] = true;
+        }
+      }
+      return;
+    }
+    for(int idx = 0; idx < nFanins; idx++) {
+      lit x = pBdd->Const1();
+      IncRef(x);
+      for(unsigned idx2 = idx + 1; idx2 < nFanins; idx2++) {
+        int fi = pNtk->GetFanin(id, idx2);
+        bool c = pNtk->GetCompl(id, idx2);
+        Assign(x, pBdd->And(x, pBdd->LitNotCond(vFs[fi], c)));
+      }
+      Assign(x, pBdd->Or(pBdd->LitNot(x), vGs[id]));
+      if(!pBdd->LitIsEq(vvCs[id][idx], x)) {
+        Assign(vvCs[id][idx], x);
+        int fi = pNtk->GetFanin(id, idx);
+        vGUpdates[fi] = true;
+      }
+      DecRef(x);
+    }
+  }
+
+  template <typename Ntk>
+  void BddAnalyzer<Ntk>::CspfNode(int id) {
+    if(!vGUpdates[id] && !vCUpdates[id]) {
+      return;
+    }
+    if(vGUpdates[id]) {
+      if(nVerbose) {
+        std::cout << "computing node " << id << " G " << std::endl;
+      }
+      if(ComputeG(id)) {
+        vCUpdates[id] = true;
+      }
+      vGUpdates[id] = false;
+    }
+    if(vCUpdates[id]) {
+      if(nVerbose) {
+        std::cout << "computing node " << id << " C " << std::endl;
+      }
+      ComputeC(id);
+      vCUpdates[id] = false;
+    }
+  }
+  
+  template <typename Ntk>
+  void BddAnalyzer<Ntk>::Cspf(int id, bool fC) {
+    if(id != -1) {
+      pNtk->ForEachTfoReverse(id, false, [&](int fo) {
+        CspfNode(fo);
+      });
+      bool fCUpdate = vCUpdates[id];
+      if(!fC) {
+        vCUpdates[id] = false;
+      }
+      CspfNode(id);
+      if(!fC) {
+        vCUpdates[id] = fCUpdate;
+      }
+    } else {
+      pNtk->ForEachIntReverse([&](int id) {
+        CspfNode(id);
+      });
+    }
+  }
+
+  /* }}} */
+
+  /* {{{ Save & load */
+
+  template <typename Ntk>
+  void BddAnalyzer<Ntk>::Save(int slot) {
+    if(slot >= vBackups.size()) {
+      vBackups.resize(slot + 1);
+    }
+    vBackups[slot].target = target;
+    CopyVec(vBackups[slot].vFs, vFs);
+    CopyVec(vBackups[slot].vGs, vGs);
+    CopyVecVec(vBackups[slot].vvCs, vvCs);
+    vBackups[slot].vUpdates = vUpdates;
+    vBackups[slot].vGUpdates = vGUpdates;
+    vBackups[slot].vCUpdates = vCUpdates;
+  }
+
+  template <typename Ntk>
+  void BddAnalyzer<Ntk>::Load(int slot) {
+    assert(slot < vBackups.size());
+    target = vBackups[slot].target;
+    CopyVec(vFs, vBackups[slot].vFs);
+    CopyVec(vGs, vBackups[slot].vGs);
+    CopyVecVec(vvCs, vBackups[slot].vvCs);
+    vUpdates = vBackups[slot].vUpdates;
+    vGUpdates = vBackups[slot].vGUpdates;
+    vCUpdates = vBackups[slot].vCUpdates;
+  }
+
+  template <typename Ntk>
+  void BddAnalyzer<Ntk>::PopBack() {
+    assert(!vBackups.empty());
+    DelVec(vBackups.back().vFs);
+    DelVec(vBackups.back().vGs);
+    DelVecVec(vBackups.back().vvCs);
+    vBackups.pop_back();
+  }
+
+  /* }}} Save & load end */
+  
+  /* {{{ Constructors */
+
+  template <typename Ntk>
+  BddAnalyzer<Ntk>::BddAnalyzer() :
+    pNtk(NULL),
+    nVerbose(0),
+    pBdd(NULL),
+    target(-1),
+    fResim(false) {
+  }
+  
+  template <typename Ntk>
+  BddAnalyzer<Ntk>::BddAnalyzer(Ntk *pNtk, Parameter const *pPar) :
+    pNtk(pNtk),
+    nVerbose(pPar->nAnalyzerVerbose),
+    target(-1),
+    fResim(false) {
+    NewBdd::Param Par;
+    Par.nObjsMaxLog = 25;
+    Par.nCacheMaxLog = 20;
+    Par.fCountOnes = true;
+    Par.nGbc = 1;
+    Par.nReo = 4000;
+    pBdd = new NewBdd::Man(pNtk->GetNumPis(), Par);
+    Allocate();
+    Assign(vFs[0], pBdd->Const0());
+    int idx = 0;
+    pNtk->ForEachPi([&](int id) {
+      Assign(vFs[id], pBdd->IthVar(idx));
+      idx++;
+    });
+    pNtk->ForEachInt([&](int id) {
+      vUpdates[id] = true;
+    });
+    Simulate();
+    pBdd->Reorder();
+    pBdd->TurnOffReo();
+    pNtk->ForEachInt([&](int id) {
+      vvCs[id].resize(pNtk->GetNumFanins(id), LitMax);
+    });
+    pNtk->ForEachPo([&](int id) {
+      vvCs[id].resize(1, LitMax);
+      Assign(vvCs[id][0], pBdd->Const0());
+      int fi = pNtk->GetFanin(id, 0);
+      vGUpdates[fi]  = true;
+    });
+    pNtk->AddCallback(std::bind(&BddAnalyzer<Ntk>::ActionCallback, this, std::placeholders::_1));
+  }
+
+  template <typename Ntk>
+  void BddAnalyzer<Ntk>::UpdateNetwork(Ntk *pNtk_, bool fSame) {
+    // clear
+    while(!vBackups.empty()) {
+      PopBack();
+    }
+    DelVec(vFs);
+    DelVec(vGs);
+    DelVecVec(vvCs);
+    pNtk = pNtk_;
+    target = -1;
+    fResim = false;
+    vUpdates.clear();
+    vGUpdates.clear();
+    vCUpdates.clear();
+    // alloc
+    bool fUseReo = false;
+    if(pBdd->GetNumVars() != pNtk->GetNumPis()) {
+      // need to reset manager
+      delete pBdd;
+      NewBdd::Param Par;
+      Par.nObjsMaxLog = 25;
+      Par.nCacheMaxLog = 20;
+      Par.fCountOnes = true;
+      Par.nGbc = 1;
+      Par.nReo = 4000;
+      pBdd = new NewBdd::Man(pNtk->GetNumPis(), Par);
+      fUseReo = true;
+    } else if(!fSame) {
+      // turning on reordering if network function changed
+      pBdd->TurnOnReo();
+      fUseReo = true;
+    }
+    Allocate();
+    // prepare
+    Assign(vFs[0], pBdd->Const0());
+    int idx = 0;
+    pNtk->ForEachPi([&](int id) {
+      Assign(vFs[id], pBdd->IthVar(idx));
+      idx++;
+    });
+    pNtk->ForEachInt([&](int id) {
+      vUpdates[id] = true;
+    });
+    Simulate();
+    if(fUseReo) {
+      pBdd->Reorder();
+      pBdd->TurnOffReo();
+    }
+    pNtk->ForEachInt([&](int id) {
+      vvCs[id].resize(pNtk->GetNumFanins(id), LitMax);
+    });
+    pNtk->ForEachPo([&](int id) {
+      vvCs[id].resize(1, LitMax);
+      Assign(vvCs[id][0], pBdd->Const0());
+      int fi = pNtk->GetFanin(id, 0);
+      vGUpdates[fi]  = true;
+    });
+    pNtk->AddCallback(std::bind(&BddAnalyzer<Ntk>::ActionCallback, this, std::placeholders::_1));
+  }
+  
+  template <typename Ntk>
+  BddAnalyzer<Ntk>::~BddAnalyzer() {
+    while(!vBackups.empty()) {
+      PopBack();
+    }
+    DelVec(vFs);
+    DelVec(vGs);
+    DelVecVec(vvCs);
+    if(pBdd) {
+      pBdd->PrintStats();
+    }
+    delete pBdd;
+  }
+
+  /* }}} */
+
+  /* {{{ Checks */
+  
+  template <typename Ntk>
+  bool BddAnalyzer<Ntk>::CheckRedundancy(int id, int idx) {
+    if(fResim) {
+      Simulate();
+    }
+    Cspf(id);
+    switch(pNtk->GetNodeType(id)) {
+    case AND: {
+      int fi = pNtk->GetFanin(id, idx);
+      bool c = pNtk->GetCompl(id, idx);
+      lit x = pBdd->Or(pBdd->LitNotCond(vFs[fi], c), vvCs[id][idx]);
+      if(pBdd->IsConst1(x)) {
+        if(nVerbose) {
+          std::cout << "node " << id << " fanin " << (pNtk->GetCompl(id, idx)? "!": "") << pNtk->GetFanin(id, idx) << " index " << idx << " is redundant" << std::endl;
+        }
+        return true;
+      }
+      break;
+    }
+    default:
+      assert(0);
+    }
+    if(nVerbose) {
+      std::cout << "node " << id << " fanin " << (pNtk->GetCompl(id, idx)? "!": "") << pNtk->GetFanin(id, idx) << " index " << idx << " is NOT redundant" << std::endl;
+    }
+    return false;
+  }
+  
+  template <typename Ntk>
+  bool BddAnalyzer<Ntk>::CheckFeasibility(int id, int fi, bool c) {
+    if(target != id) { // simualte if there has been update in non-tfo of this node
+      Simulate();
+      target = id;
+    }
+    Cspf(id, false);
+    switch(pNtk->GetNodeType(id)) {
+    case AND: {
+      lit x = pBdd->Or(pBdd->LitNot(vFs[id]), vGs[id]);
+      IncRef(x);
+      lit y = pBdd->Or(x, pBdd->LitNotCond(vFs[fi], c));
+      DecRef(x);
+      if(pBdd->IsConst1(y)) {
+        if(nVerbose) {
+          std::cout << "node " << id << " fanin " << (c? "!": "") << fi << " is feasible" << std::endl;
+        }
+        return true;
+      }
+      break;
+    }
+    default:
+      assert(0);
+    }
+    if(nVerbose) {
+      std::cout << "node " << id << " fanin " << (c? "!": "") << fi << " is NOT feasible" << std::endl;
+    }
+    return false;
+  }
+
+  /* }}} */
+  
+}
+
+ABC_NAMESPACE_CXX_HEADER_END
diff --git a/src/opt/rrr/rrrBddMspfAnalyzer.h b/src/opt/rrr/rrrBddMspfAnalyzer.h
new file mode 100644
index 000000000..6c3b057ed
--- /dev/null
+++ b/src/opt/rrr/rrrBddMspfAnalyzer.h
@@ -0,0 +1,761 @@
+#pragma once
+
+#include <iostream>
+#include <vector>
+
+#include <aig/gia/giaNewBdd.h>
+
+#include "rrrParameter.h"
+#include "rrrTypes.h"
+
+ABC_NAMESPACE_CXX_HEADER_START
+
+namespace rrr {
+
+  template <typename Ntk>
+  class BddMspfAnalyzer {
+  private:
+    // aliases
+    using lit = int;
+    static constexpr lit LitMax = 0xffffffff;
+
+    // pointer to network
+    Ntk *pNtk;
+    
+    // parameters
+    int nVerbose;
+    
+    // data
+    NewBdd::Man *pBdd;
+    std::vector<lit> vFs;
+    std::vector<lit> vGs;
+    std::vector<std::vector<lit>> vvCs;
+    bool fUpdate;
+    std::vector<bool> vUpdates;
+    std::vector<bool> vGUpdates;
+    std::vector<bool> vCUpdates;
+    std::vector<bool> vVisits;
+    
+    // backups
+    std::vector<BddMspfAnalyzer> vBackups;
+
+    // BDD utils
+    void IncRef(lit x) const;
+    void DecRef(lit x) const;
+    void Assign(lit &x, lit y) const;
+    void CopyVec(std::vector<lit> &x, std::vector<lit> const &y) const;
+    void DelVec(std::vector<lit> &v) const;
+    void CopyVecVec(std::vector<std::vector<lit>> &x, std::vector<std::vector<lit>> const &y) const;
+    void DelVecVec(std::vector<std::vector<lit>> &v) const;
+    lit  Xor(lit x, lit y) const;
+
+    // callback
+    void ActionCallback(Action const &action);
+
+    // allocation
+    void Allocate();
+
+    // simulation
+    bool SimulateNode(int id, std::vector<lit> &v) const;
+    void Simulate();
+
+    // PF computation
+    bool ComputeReconvergentG(int id);
+    bool ComputeG(int id);
+    void ComputeC(int id);
+    bool ComputeCDebug(int id);
+    void MspfNode(int id);
+    void Mspf(int id = -1, bool fC = true);
+
+    // save & load
+    void Save(int slot);
+    void Load(int slot);
+    void PopBack();
+    
+  public:
+    // constructors
+    BddMspfAnalyzer();
+    BddMspfAnalyzer(Ntk *pNtk, Parameter const *pPar);
+    ~BddMspfAnalyzer();
+    void UpdateNetwork(Ntk *pNtk_, bool fSame);
+
+    // checks
+    bool CheckRedundancy(int id, int idx);
+    bool CheckFeasibility(int id, int fi, bool c);
+  };
+  
+  /* {{{ BDD utils */
+  
+  template <typename Ntk>
+  inline void BddMspfAnalyzer<Ntk>::IncRef(lit x) const {
+    if(x != LitMax) {
+      pBdd->IncRef(x);
+    }
+  }
+  
+  template <typename Ntk>
+  inline void BddMspfAnalyzer<Ntk>::DecRef(lit x) const {
+    if(x != LitMax) {
+      pBdd->DecRef(x);
+    }
+  }
+  
+  template <typename Ntk>
+  inline void BddMspfAnalyzer<Ntk>::Assign(lit &x, lit y) const {
+    DecRef(x);
+    x = y;
+    IncRef(x);
+  }
+
+  template <typename Ntk>
+  inline void BddMspfAnalyzer<Ntk>::CopyVec(std::vector<lit> &x, std::vector<lit> const &y) const {
+    for(size_t i = y.size(); i < x.size(); i++) {
+      DecRef(x[i]);
+    }
+    x.resize(y.size(), LitMax);
+    for(size_t i = 0; i < y.size(); i++) {
+      if(x[i] != y[i]) {
+        DecRef(x[i]);
+        x[i] = y[i];
+        IncRef(x[i]);
+      }
+    }
+  }
+
+  template <typename Ntk>
+  inline void BddMspfAnalyzer<Ntk>::DelVec(std::vector<lit> &v) const {
+    for(lit x: v) {
+      DecRef(x);
+    }
+    v.clear();
+  }
+
+  template <typename Ntk>
+  inline void BddMspfAnalyzer<Ntk>::CopyVecVec(std::vector<std::vector<lit>> &x, std::vector<std::vector<lit>> const &y) const {
+    for(size_t i = y.size(); i < x.size(); i++) {
+      DelVec(x[i]);
+    }
+    x.resize(y.size());
+    for(size_t i = 0; i < y.size(); i++) {
+      CopyVec(x[i], y[i]);
+    }
+  }
+
+  template <typename Ntk>
+  inline void BddMspfAnalyzer<Ntk>::DelVecVec(std::vector<std::vector<lit>> &v) const {
+    for(size_t i = 0; i < v.size(); i++) {
+      DelVec(v[i]);
+    }
+    v.clear();
+  }
+  
+  template <typename Ntk>
+  inline typename BddMspfAnalyzer<Ntk>::lit BddMspfAnalyzer<Ntk>::Xor(lit x, lit y) const {
+    lit f = pBdd->And(x, pBdd->LitNot(y));
+    pBdd->IncRef(f);
+    lit g = pBdd->And(pBdd->LitNot(x), y);
+    pBdd->IncRef(g);
+    lit r = pBdd->Or(f, g);
+    pBdd->DecRef(f);
+    pBdd->DecRef(g);
+    return r;
+  }
+  
+  /* }}} */
+
+  /* {{{ Callback */
+
+  template <typename Ntk>
+  void BddMspfAnalyzer<Ntk>::ActionCallback(Action const &action) {
+    switch(action.type) {
+    case REMOVE_FANIN:
+      fUpdate = true;
+      std::fill(vVisits.begin(), vVisits.end(), false);
+      vUpdates[action.id] = true;
+      vCUpdates[action.id] = true;
+      vGUpdates[action.fi] = true;
+      DecRef(vvCs[action.id][action.idx]);
+      vvCs[action.id].erase(vvCs[action.id].begin() + action.idx);
+      break;
+    case REMOVE_UNUSED:
+      if(vGUpdates[action.id] || vCUpdates[action.id]) {
+        for(int fi: action.vFanins) {
+          vGUpdates[fi] = true;
+        }
+      }
+      Assign(vGs[action.id], LitMax);
+      DelVec(vvCs[action.id]);
+      break;
+    case REMOVE_BUFFER:
+      if(vUpdates[action.id]) {
+        fUpdate = true;
+        for(int fo: action.vFanouts) {
+          vUpdates[fo] = true;
+          vCUpdates[fo] = true;
+        }
+      }
+      if(vGUpdates[action.id] || vCUpdates[action.id]) {
+        vGUpdates[action.fi] = true;
+      }
+      Assign(vGs[action.id], LitMax);
+      DelVec(vvCs[action.id]);
+      break;
+    case REMOVE_CONST:
+      if(vUpdates[action.id]) {
+        fUpdate = true;
+        for(int fo: action.vFanouts) {
+          vUpdates[fo] = true;
+          vCUpdates[fo] = true;
+        }
+      }
+      for(int fi: action.vFanins) {
+        vGUpdates[fi] = true;
+      }
+      break;
+    case ADD_FANIN:
+      fUpdate = true;
+      std::fill(vVisits.begin(), vVisits.end(), false);
+      vUpdates[action.id] = true;
+      vCUpdates[action.id] = true;
+      vvCs[action.id].insert(vvCs[action.id].begin() + action.idx, LitMax);
+      break;
+    case TRIVIAL_COLLAPSE: {
+      if(vGUpdates[action.fi] || vCUpdates[action.fi]) {
+        vCUpdates[action.id] = true;
+      }
+      std::vector<lit>::iterator it = vvCs[action.id].begin() + action.idx;
+      DecRef(*it);
+      it = vvCs[action.id].erase(it);
+      vvCs[action.id].insert(it,  vvCs[action.fi].begin(), vvCs[action.fi].end());
+      vvCs[action.fi].clear();
+      Assign(vFs[action.fi], LitMax);
+      Assign(vGs[action.fi], LitMax);
+      break;
+    }
+    case TRIVIAL_DECOMPOSE: {
+      Allocate();
+      SimulateNode(action.fi, vFs);
+      assert(vGs[action.fi] == LitMax);
+      std::vector<lit>::iterator it = vvCs[action.id].begin() + action.idx;
+      assert(vvCs[action.fi].empty());
+      vvCs[action.fi].insert(vvCs[action.fi].begin(), it, vvCs[action.id].end());
+      vvCs[action.id].erase(it, vvCs[action.id].end());
+      assert(vvCs[action.id].size() == action.idx);
+      if(!vGUpdates[action.id] && !vCUpdates[action.id]) {
+        // recompute here only when updates are unlikely to happen
+        if(pBdd->IsConst1(vGs[action.id])) {
+          Assign(vGs[action.fi], pBdd->Const1());
+        } else {
+          lit x = pBdd->Const1();
+          IncRef(x);
+          for(int idx2 = 0; idx2 < action.idx; idx2++) {
+            int fi = pNtk->GetFanin(action.id, idx2);
+            bool c = pNtk->GetCompl(action.id, idx2);
+            Assign(x, pBdd->And(x, pBdd->LitNotCond(vFs[fi], c)));
+          }
+          Assign(vGs[action.fi], pBdd->Or(pBdd->LitNot(x), vGs[action.id]));
+          DecRef(x);
+        }
+      }
+      vvCs[action.id].resize(action.idx + 1, LitMax);
+      Assign(vvCs[action.id][action.idx], vGs[action.fi]);
+      vUpdates[action.fi] = false;
+      vGUpdates[action.fi] = false;
+      vCUpdates[action.fi] = false;
+      break;
+    }
+    case SORT_FANINS: {
+      std::vector<lit> vCs = vvCs[action.id];
+      vvCs[action.id].clear();
+      for(int index: action.vIndices) {
+        vvCs[action.id].push_back(vCs[index]);
+      }
+      break;
+    }
+    case SAVE:
+      Save(action.idx);
+      break;
+    case LOAD:
+      Load(action.idx);
+      break;
+    case POP_BACK:
+      PopBack();
+      break;
+    default:
+      assert(0);
+    }
+  }
+  
+  /* }}} */
+
+  /* {{{ Allocation */
+
+  template <typename Ntk>
+  void BddMspfAnalyzer<Ntk>::Allocate() {
+    int nNodes = pNtk->GetNumNodes();
+    vFs.resize(nNodes, LitMax);
+    vGs.resize(nNodes, LitMax);
+    vvCs.resize(nNodes);
+    vUpdates.resize(nNodes);
+    vGUpdates.resize(nNodes);
+    vCUpdates.resize(nNodes);
+    vVisits.resize(nNodes);
+  }
+
+  /* }}} */
+  
+  /* {{{ Simulation */
+  
+  template <typename Ntk>
+  inline bool BddMspfAnalyzer<Ntk>::SimulateNode(int id, std::vector<lit> &v) const {
+    if(nVerbose) {
+      std::cout << "simulating node " << id << std::endl;
+    }
+    lit x = pBdd->Const1();
+    IncRef(x);
+    pNtk->ForEachFanin(id, [&](int fi, bool c) {
+      Assign(x, pBdd->And(x, pBdd->LitNotCond(v[fi], c)));
+    });
+    if(pBdd->LitIsEq(x, v[id])) {
+      DecRef(x);
+      return false;
+    }
+    Assign(v[id], x);
+    DecRef(x);
+    return true;
+  }
+  
+  template <typename Ntk>
+  void BddMspfAnalyzer<Ntk>::Simulate() {
+    if(nVerbose) {
+      std::cout << "symbolic simulation with BDD" << std::endl;
+    }
+    pNtk->ForEachInt([&](int id) {
+      if(vUpdates[id]) {
+        if(SimulateNode(id, vFs)) {
+          pNtk->ForEachFanout(id, false, [&](int fo, bool c) {
+            vUpdates[fo] = true;
+            vCUpdates[fo] = true;
+          });
+        }
+        vUpdates[id] = false;
+      }
+    });
+  }
+
+  /* }}} */
+
+  /* {{{ MSPF computation */
+
+  template <typename pNtk>
+  inline bool BddMspfAnalyzer<pNtk>::ComputeReconvergentG(int id) {
+    std::vector<lit> v;
+    CopyVec(v, vFs);
+    v[id] = pBdd->LitNot(v[id]);
+    pNtk->ForEachTfoUpdate(id, false, [&](int fo) {
+      return SimulateNode(fo, v);
+    });
+    lit x = pBdd->Const1();
+    IncRef(x);
+    pNtk->ForEachPoDriver([&](int fi, bool c) {
+      lit y = Xor(vFs[fi], v[fi]);
+      IncRef(y);
+      Assign(x, pBdd->And(x, pBdd->LitNot(y)));
+      DecRef(y);
+    });
+    if(pBdd->LitIsEq(vGs[id], x)) {
+      DecRef(x);
+      DelVec(v);
+      return false;
+    }
+    Assign(vGs[id], x);
+    DecRef(x);
+    DelVec(v);
+    return true;
+  }
+  
+  template <typename pNtk>
+  inline bool BddMspfAnalyzer<pNtk>::ComputeG(int id) {
+    if(pNtk->GetNumFanouts(id) == 0) {
+      if(pBdd->IsConst1(vGs[id])) {
+        return false;
+      }
+      Assign(vGs[id], pBdd->Const1());
+      return true;
+    }
+    lit x = pBdd->Const1();
+    IncRef(x);
+    pNtk->ForEachFanoutRidx(id, true, [&](int fo, bool c, int idx) {
+      Assign(x, pBdd->And(x, vvCs[fo][idx]));
+    });
+    if(pBdd->LitIsEq(vGs[id], x)) {
+      DecRef(x);
+      return false;
+    }
+    Assign(vGs[id], x);
+    DecRef(x);
+    return true;
+  }
+
+  template <typename pNtk>
+  inline void BddMspfAnalyzer<pNtk>::ComputeC(int id) {
+    int nFanins = pNtk->GetNumFanins(id);
+    assert(vvCs[id].size() == nFanins);
+    if(pBdd->IsConst1(vGs[id])) {
+      for(int idx = 0; idx < nFanins; idx++) {
+        if(!pBdd->IsConst1(vvCs[id][idx])) {
+          Assign(vvCs[id][idx], pBdd->Const1());
+          int fi = pNtk->GetFanin(id, idx);
+          vGUpdates[fi] = true;
+        }
+      }
+      return;
+    }
+    for(int idx = 0; idx < nFanins; idx++) {
+      lit x = pBdd->Const1();
+      IncRef(x);
+      for(int idx2 = 0; idx2 < nFanins; idx2++) {
+        if(idx2 != idx) {
+          int fi = pNtk->GetFanin(id, idx2);
+          bool c = pNtk->GetCompl(id, idx2);
+          Assign(x, pBdd->And(x, pBdd->LitNotCond(vFs[fi], c)));
+        }
+      }
+      Assign(x, pBdd->Or(pBdd->LitNot(x), vGs[id]));
+      if(!pBdd->LitIsEq(vvCs[id][idx], x)) {
+        Assign(vvCs[id][idx], x);
+        int fi = pNtk->GetFanin(id, idx);
+        vGUpdates[fi] = true;
+      }
+      DecRef(x);
+    }
+  }
+  
+  template <typename pNtk>
+  inline bool BddMspfAnalyzer<pNtk>::ComputeCDebug(int id) {
+    bool fUpdated = false;
+    int nFanins = pNtk->GetNumFanins(id);
+    assert(vvCs[id].size() == nFanins);
+    if(pBdd->IsConst1(vGs[id])) {
+      for(int idx = 0; idx < nFanins; idx++) {
+        if(!pBdd->IsConst1(vvCs[id][idx])) {
+          Assign(vvCs[id][idx], pBdd->Const1());
+          int fi = pNtk->GetFanin(id, idx);
+          vGUpdates[fi] = true;
+          fUpdated = true;
+        }
+      }
+      return fUpdated;
+    }
+    for(int idx = 0; idx < nFanins; idx++) {
+      lit x = pBdd->Const1();
+      IncRef(x);
+      for(int idx2 = 0; idx2 < nFanins; idx2++) {
+        if(idx2 != idx) {
+          int fi = pNtk->GetFanin(id, idx2);
+          bool c = pNtk->GetCompl(id, idx2);
+          Assign(x, pBdd->And(x, pBdd->LitNotCond(vFs[fi], c)));
+        }
+      }
+      Assign(x, pBdd->Or(pBdd->LitNot(x), vGs[id]));
+      if(!pBdd->LitIsEq(vvCs[id][idx], x)) {
+        Assign(vvCs[id][idx], x);
+        int fi = pNtk->GetFanin(id, idx);
+        vGUpdates[fi] = true;
+        fUpdated = true;
+      }
+      DecRef(x);
+    }
+    return fUpdated;
+  }
+
+  template <typename pNtk>
+  void BddMspfAnalyzer<pNtk>::MspfNode(int id) {
+    if(vGUpdates[id]) {
+      if(pNtk->IsReconvergent(id)) {
+        if(nVerbose) {
+          std::cout << "computing reconvergent node " << id << " G " << std::endl;
+        }
+        if(ComputeReconvergentG(id)) {
+          vCUpdates[id] = true;
+        }
+      } else {
+        if(nVerbose) {
+          std::cout << "computing node " << id << " G " << std::endl;
+        }
+        if(ComputeG(id)) {
+          vCUpdates[id] = true;
+        }
+      }
+      vGUpdates[id] = false;
+    } else if(!vVisits[id] && pNtk->IsReconvergent(id)) {
+      if(nVerbose) {
+        std::cout << "computing unvisited reconvergent node " << id << " G " << std::endl;
+      }
+      if(ComputeReconvergentG(id)) {
+        vCUpdates[id] = true;
+      }
+    }
+    if(vCUpdates[id]) {
+      if(nVerbose) {
+        std::cout << "computing node " << id << " C " << std::endl;
+      }
+      ComputeC(id);
+      vCUpdates[id] = false;
+    }
+    vVisits[id] = true;
+  }
+
+  template <typename pNtk>
+  void BddMspfAnalyzer<pNtk>::Mspf(int id, bool fC) {
+    if(id != -1) {
+      pNtk->ForEachTfoReverse(id, false, [&](int fo) {
+        MspfNode(fo);
+      });
+      bool fCUpdate = vCUpdates[id];
+      if(!fC) {
+        vCUpdates[id] = false;
+      }
+      MspfNode(id);
+      if(!fC) {
+        vCUpdates[id] = fCUpdate;
+      }
+    } else {
+      pNtk->ForEachIntReverse([&](int id) {
+        MspfNode(id);
+      });
+    }
+  }
+
+  /* }}} */
+
+  /* {{{ Save & load */
+
+  template <typename Ntk>
+  void BddMspfAnalyzer<Ntk>::Save(int slot) {
+    if(slot >= vBackups.size()) {
+      vBackups.resize(slot + 1);
+    }
+    CopyVec(vBackups[slot].vFs, vFs);
+    CopyVec(vBackups[slot].vGs, vGs);
+    CopyVecVec(vBackups[slot].vvCs, vvCs);
+    vBackups[slot].fUpdate = fUpdate;
+    vBackups[slot].vUpdates = vUpdates;
+    vBackups[slot].vGUpdates = vGUpdates;
+    vBackups[slot].vCUpdates = vCUpdates;
+    vBackups[slot].vVisits = vVisits;
+  }
+
+  template <typename Ntk>
+  void BddMspfAnalyzer<Ntk>::Load(int slot) {
+    assert(slot < vBackups.size());
+    CopyVec(vFs, vBackups[slot].vFs);
+    CopyVec(vGs, vBackups[slot].vGs);
+    CopyVecVec(vvCs, vBackups[slot].vvCs);
+    fUpdate = vBackups[slot].fUpdate;
+    vUpdates = vBackups[slot].vUpdates;
+    vGUpdates = vBackups[slot].vGUpdates;
+    vCUpdates = vBackups[slot].vCUpdates;
+    vVisits = vBackups[slot].vVisits;
+  }
+
+  template <typename Ntk>
+  void BddMspfAnalyzer<Ntk>::PopBack() {
+    assert(!vBackups.empty());
+    DelVec(vBackups.back().vFs);
+    DelVec(vBackups.back().vGs);
+    DelVecVec(vBackups.back().vvCs);
+    vBackups.pop_back();
+  }
+
+  /* }}} */
+  
+  /* {{{ Constructor */
+
+  template <typename Ntk>
+  BddMspfAnalyzer<Ntk>::BddMspfAnalyzer() :
+    pNtk(NULL),
+    nVerbose(0),
+    pBdd(NULL) {
+  }
+  
+  template <typename Ntk>
+  BddMspfAnalyzer<Ntk>::BddMspfAnalyzer(Ntk *pNtk, Parameter const *pPar) :
+    pNtk(pNtk),
+    nVerbose(pPar->nAnalyzerVerbose),
+    fUpdate(false) {
+    NewBdd::Param Par;
+    Par.nObjsMaxLog = 25;
+    Par.nCacheMaxLog = 20;
+    Par.fCountOnes = true;
+    Par.nGbc = 1;
+    Par.nReo = 4000;
+    pBdd = new NewBdd::Man(pNtk->GetNumPis(), Par);
+    Allocate();
+    Assign(vFs[0], pBdd->Const0());
+    int idx = 0;
+    pNtk->ForEachPi([&](int id) {
+      Assign(vFs[id], pBdd->IthVar(idx));
+      idx++;
+    });
+    pNtk->ForEachInt([&](int id) {
+      vUpdates[id] = true;
+    });
+    Simulate();
+    pBdd->Reorder();
+    pBdd->TurnOffReo();
+    pNtk->ForEachInt([&](int id) {
+      vvCs[id].resize(pNtk->GetNumFanins(id), LitMax);
+    });
+    pNtk->ForEachPo([&](int id) {
+      vvCs[id].resize(1, LitMax);
+      Assign(vvCs[id][0], pBdd->Const0());
+      int fi = pNtk->GetFanin(id, 0);
+      vGUpdates[fi]  = true;
+    });
+    pNtk->AddCallback(std::bind(&BddMspfAnalyzer<Ntk>::ActionCallback, this, std::placeholders::_1));
+  }
+
+  template <typename Ntk>
+  void BddMspfAnalyzer<Ntk>::UpdateNetwork(Ntk *pNtk_, bool fSame) {
+    // clear
+    while(!vBackups.empty()) {
+      PopBack();
+    }
+    DelVec(vFs);
+    DelVec(vGs);
+    DelVecVec(vvCs);
+    pNtk = pNtk_;
+    fUpdate = false;
+    vUpdates.clear();
+    vGUpdates.clear();
+    vCUpdates.clear();
+    vVisits.clear();
+    // alloc
+    bool fUseReo = false;
+    if(pBdd->GetNumVars() != pNtk->GetNumPis()) {
+      // need to reset manager
+      delete pBdd;
+      NewBdd::Param Par;
+      Par.nObjsMaxLog = 25;
+      Par.nCacheMaxLog = 20;
+      Par.fCountOnes = true;
+      Par.nGbc = 1;
+      Par.nReo = 4000;
+      pBdd = new NewBdd::Man(pNtk->GetNumPis(), Par);
+      fUseReo = true;
+    } else if(!fSame) {
+      // turning on reordering if network function changed
+      pBdd->TurnOnReo();
+      fUseReo = true;
+    }
+    Allocate();
+    // prepare
+    Assign(vFs[0], pBdd->Const0());
+    int idx = 0;
+    pNtk->ForEachPi([&](int id) {
+      Assign(vFs[id], pBdd->IthVar(idx));
+      idx++;
+    });
+    pNtk->ForEachInt([&](int id) {
+      vUpdates[id] = true;
+    });
+    Simulate();
+    if(fUseReo) {
+      pBdd->Reorder();
+      pBdd->TurnOffReo();
+    }
+    pNtk->ForEachInt([&](int id) {
+      vvCs[id].resize(pNtk->GetNumFanins(id), LitMax);
+    });
+    pNtk->ForEachPo([&](int id) {
+      vvCs[id].resize(1, LitMax);
+      Assign(vvCs[id][0], pBdd->Const0());
+      int fi = pNtk->GetFanin(id, 0);
+      vGUpdates[fi]  = true;
+    });
+    pNtk->AddCallback(std::bind(&BddMspfAnalyzer<Ntk>::ActionCallback, this, std::placeholders::_1));
+  }
+  
+  template <typename Ntk>
+  BddMspfAnalyzer<Ntk>::~BddMspfAnalyzer() {
+    while(!vBackups.empty()) {
+      PopBack();
+    }
+    DelVec(vFs);
+    DelVec(vGs);
+    DelVecVec(vvCs);
+    if(pBdd) {
+      pBdd->PrintStats();
+    }
+    delete pBdd;
+  }
+
+  /* }}} */
+
+  /* {{{ Checks */
+  
+  template <typename Ntk>
+  bool BddMspfAnalyzer<Ntk>::CheckRedundancy(int id, int idx) {
+    if(fUpdate) {
+      Simulate();
+      fUpdate = false;
+    }
+    Mspf(id);
+    switch(pNtk->GetNodeType(id)) {
+    case AND: {
+      int fi = pNtk->GetFanin(id, idx);
+      bool c = pNtk->GetCompl(id, idx);
+      lit x = pBdd->Or(pBdd->LitNotCond(vFs[fi], c), vvCs[id][idx]);
+      if(pBdd->IsConst1(x)) {
+        if(nVerbose) {
+          std::cout << "node " << id << " fanin " << (pNtk->GetCompl(id, idx)? "!": "") << pNtk->GetFanin(id, idx) << " index " << idx << " is redundant" << std::endl;
+        }
+        return true;
+      }
+      break;
+    }
+    default:
+      assert(0);
+    }
+    if(nVerbose) {
+      std::cout << "node " << id << " fanin " << (pNtk->GetCompl(id, idx)? "!": "") << pNtk->GetFanin(id, idx) << " index " << idx << " is NOT redundant" << std::endl;
+    }
+    return false;
+  }
+  
+  template <typename Ntk>
+  bool BddMspfAnalyzer<Ntk>::CheckFeasibility(int id, int fi, bool c) {
+    if(fUpdate) {
+      Simulate();
+      fUpdate = false;
+    }
+    Mspf(id, false);
+    switch(pNtk->GetNodeType(id)) {
+    case AND: {
+      lit x = pBdd->Or(pBdd->LitNot(vFs[id]), vGs[id]);
+      IncRef(x);
+      lit y = pBdd->Or(x, pBdd->LitNotCond(vFs[fi], c));
+      DecRef(x);
+      if(pBdd->IsConst1(y)) {
+        if(nVerbose) {
+          std::cout << "node " << id << " fanin " << (c? "!": "") << fi << " is feasible" << std::endl;
+        }
+        return true;
+      }
+      break;
+    }
+    default:
+      assert(0);
+    }
+    if(nVerbose) {
+      std::cout << "node " << id << " fanin " << (c? "!": "") << fi << " is NOT feasible" << std::endl;
+    }
+    return false;
+  }
+
+  /* }}} */
+  
+}
+
+ABC_NAMESPACE_CXX_HEADER_END
diff --git a/src/opt/rrr/rrrOptimizer.h b/src/opt/rrr/rrrOptimizer.h
new file mode 100644
index 000000000..c114bfe3b
--- /dev/null
+++ b/src/opt/rrr/rrrOptimizer.h
@@ -0,0 +1,869 @@
+#pragma once
+
+#include <vector>
+#include <map>
+#include <iterator>
+#include <random>
+#include <numeric>
+#include <limits>
+
+#include "rrrParameter.h"
+#include "rrrTypes.h"
+
+ABC_NAMESPACE_CXX_HEADER_START
+
+namespace rrr {
+
+  template <typename Ntk, typename Ana>
+  class Optimizer {
+  private:
+    // aliases
+    using itr = std::vector<int>::iterator;
+    using citr = std::vector<int>::const_iterator;
+    using ritr = std::vector<int>::reverse_iterator;
+    using critr = std::vector<int>::const_reverse_iterator;
+
+    // pointer to network
+    Ntk *pNtk;
+    
+    // parameters
+    int nVerbose;
+    std::function<double(Ntk *)> CostFunction;
+    int nSortType;
+    int nFlow;
+    seconds nTimeout; // assigned upon Run
+
+    // data
+    Ana *pAna;
+    std::mt19937 rng;
+    std::vector<int> vTmp;
+    std::map<int, std::set<int>> mapNewFanins;
+    time_point start;
+
+    // fanin sorting data
+    std::vector<int> vRandPiOrder;
+    std::vector<double> vRandCosts;
+
+    // marks
+    int target;
+    std::vector<bool> vMarks;
+    
+    // callback
+    void ActionCallback(Action const &action);
+
+    // topology
+    void MarkTfo(int id);
+
+    // time
+    bool Timeout();
+
+    // sort fanins
+    void SetRandPiOrder();
+    void SetRandCosts();
+    void SortFanins(int id);
+
+    // reduce fanin
+    bool ReduceFanin(int id, bool fRemoveUnused = false);
+    bool ReduceFaninOneRandom(int id, bool fRemoveUnused = false);
+
+    // reduce
+    void Reduce();
+    void ReduceRandom();
+
+    // remove redundancy
+    void RemoveRedundancy();
+    void RemoveRedundancyRandom();
+
+    // addition
+    template <typename T>
+    T SingleAdd(int id, T begin, T end);
+    int  MultiAdd(int id, std::vector<int> const &vCands, int nMax = 0);
+
+    // resub
+    void SingleResub(int id, std::vector<int> const &vCands, bool fGreedy = true);
+    void MultiResub(int id, std::vector<int> const &vCands, bool fGreedy = true, int nMax = 0);
+    
+    // apply
+    void ApplyReverseTopologically(std::function<void(int)> const &func);
+    
+  public:
+    // constructors
+    Optimizer(Ntk *pNtk, Parameter const *pPar, std::function<double(Ntk *)> CostFunction);
+    ~Optimizer();
+    void UpdateNetwork(Ntk *pNtk_, bool fSame);
+
+    // run
+    void Run(seconds nTimeout_ = 0);
+    void Randomize();
+    
+  };
+
+  /* {{{ Callback */
+  
+  template <typename Ntk, typename Ana>
+  void Optimizer<Ntk, Ana>::ActionCallback(Action const &action) {
+    if(nVerbose) {
+      PrintAction(action);
+    }
+    switch(action.type) {
+    case REMOVE_FANIN:
+      if(action.id != target) {
+        target = -1;
+      }
+      break;
+    case REMOVE_UNUSED:
+      break;
+    case REMOVE_BUFFER:
+    case REMOVE_CONST:
+      if(action.id == target) {
+        target = -1;
+      }
+      break;
+    case ADD_FANIN:
+      if(action.id != target) {
+        target = -1;
+      }
+      break;
+    case TRIVIAL_COLLAPSE:
+      break;
+    case TRIVIAL_DECOMPOSE:
+      target = -1;
+      break;
+    case SORT_FANINS:
+      break;
+    case SAVE:
+      break;
+    case LOAD:
+      //target = -1; // this is not always needed, so do it manually
+      break;
+    case POP_BACK:
+      break;
+    default:
+      assert(0);
+    }
+  }
+
+  /* }}} */
+
+  /* {{{ Topology */
+  
+  template <typename Ntk, typename Ana>
+  inline void Optimizer<Ntk, Ana>::MarkTfo(int id) {
+    // includes id itself
+    if(id == target) {
+      return;
+    }
+    target = id;
+    vMarks.clear();
+    vMarks.resize(pNtk->GetNumNodes());
+    vMarks[id] = true;
+    pNtk->ForEachTfo(id, false, [&](int fo) {
+      vMarks[fo] = true;
+    });
+  }
+
+  /* }}} */
+
+  /* {{{ Time */
+  
+  template <typename Ntk, typename Ana>
+  inline bool Optimizer<Ntk, Ana>::Timeout() {
+    if(nTimeout) {
+      time_point current = GetCurrentTime();
+      if(DurationInSeconds(start, current) > nTimeout) {
+        return true;
+      }
+    }
+    return false;
+  }
+
+  /* }}} */
+
+  /* {{{ Sort fanins */
+
+  template <typename Ntk, typename Ana>
+  inline void Optimizer<Ntk, Ana>::SetRandPiOrder() {
+    assert(vRandPiOrder.size() <= (std::vector<int>::size_type)std::numeric_limits<int>::max());
+    if((int)vRandPiOrder.size() != pNtk->GetNumPis()) {
+      vRandPiOrder.clear();
+      vRandPiOrder.resize(pNtk->GetNumPis());
+      std::iota(vRandPiOrder.begin(), vRandPiOrder.end(), 0);
+      std::shuffle(vRandPiOrder.begin(), vRandPiOrder.end(), rng);      
+    }    
+  }
+
+  template <typename Ntk, typename Ana>
+  inline void Optimizer<Ntk, Ana>::SetRandCosts() {
+    std::uniform_real_distribution<> dis(std::numeric_limits<double>::lowest(), std::numeric_limits<double>::max());
+    assert(vRandCosts.size() <= (std::vector<int>::size_type)std::numeric_limits<int>::max());
+    while((int)vRandCosts.size() < pNtk->GetNumNodes()) {
+      vRandCosts.push_back(dis(rng));
+    }
+  }
+
+  template <typename Ntk, typename Ana>
+  inline void Optimizer<Ntk, Ana>::SortFanins(int id) {
+    switch(nSortType) {
+    case 0: // no sorting
+      break;
+    case 1: // prioritize internals
+      pNtk->SortFanins(id, [&](int i, bool ci, int j, bool cj) {
+        return !pNtk->IsPi(i) && pNtk->IsPi(j);
+      });
+      break;
+    case 2: // prioritize internals with (reversely) sorted PIs
+      pNtk->SortFanins(id, [&](int i, bool ci, int j, bool cj) {
+        if(pNtk->IsPi(i) && pNtk->IsPi(j)) {
+          return pNtk->GetPiIndex(i) > pNtk->GetPiIndex(j);
+        }
+        return pNtk->IsPi(j);
+      });
+      break;
+    case 3: // prioritize internals with random PI order
+      SetRandPiOrder();
+      pNtk->SortFanins(id, [&](int i, bool ci, int j, bool cj) {
+        if(pNtk->IsPi(i) && pNtk->IsPi(j)) {
+          return vRandPiOrder[pNtk->GetPiIndex(i)] > vRandPiOrder[pNtk->GetPiIndex(j)];
+        }
+        return pNtk->IsPi(j);
+      });
+      break;
+    case 4: // smaller fanout takes larger cost
+      pNtk->SortFanins(id, [&](int i, bool ci, int j, bool cj) {
+        return pNtk->GetNumFanouts(i) < pNtk->GetNumFanouts(j);
+      });
+      break;
+    case 5: // fanout + PI
+      pNtk->SortFanins(id, [&](int i, bool ci, int j, bool cj) {
+        if(pNtk->IsPi(i) && !pNtk->IsPi(j)) {
+          return false;
+        }
+        if(!pNtk->IsPi(i) && pNtk->IsPi(j)) {
+          return true;
+        }
+        return pNtk->GetNumFanouts(i) < pNtk->GetNumFanouts(j);
+      });
+      break;
+    case 6: // fanout + sorted PI
+      pNtk->SortFanins(id, [&](int i, bool ci, int j, bool cj) {
+        if(pNtk->IsPi(i) && pNtk->IsPi(j)) {
+          return pNtk->GetPiIndex(i) > pNtk->GetPiIndex(j);
+        }
+        if(pNtk->IsPi(i)) {
+          return false;
+        }
+        if(pNtk->IsPi(j)) {
+          return true;
+        }
+        return pNtk->GetNumFanouts(i) < pNtk->GetNumFanouts(j);
+      });
+      break;
+    case 7: // fanout + random PI
+      SetRandPiOrder();
+      pNtk->SortFanins(id, [&](int i, bool ci, int j, bool cj) {
+        if(pNtk->IsPi(i) && pNtk->IsPi(j)) {
+          return vRandPiOrder[pNtk->GetPiIndex(i)] > vRandPiOrder[pNtk->GetPiIndex(j)];
+        }
+        if(pNtk->IsPi(i)) {
+          return false;
+        }
+        if(pNtk->IsPi(j)) {
+          return true;
+        }
+        return pNtk->GetNumFanouts(i) < pNtk->GetNumFanouts(j);
+      });
+      break;
+    case 8: // reverse topological order + PI
+      pNtk->SortFanins(id, [&](int i, bool ci, int j, bool cj) {
+        if(pNtk->IsPi(i) && pNtk->IsPi(j)) {
+          return false;
+        }
+        if(pNtk->IsPi(i)) {
+          return false;
+        }
+        if(pNtk->IsPi(j)) {
+          return true;
+        }
+        return pNtk->GetIntIndex(i) > pNtk->GetIntIndex(j);
+      });
+      break;
+    case 9: // reverse topological order + sorted PI
+      pNtk->SortFanins(id, [&](int i, bool ci, int j, bool cj) {
+        if(pNtk->IsPi(i) && pNtk->IsPi(j)) {
+          return pNtk->GetPiIndex(i) > pNtk->GetPiIndex(j);
+        }
+        if(pNtk->IsPi(i)) {
+          return false;
+        }
+        if(pNtk->IsPi(j)) {
+          return true;
+        }
+        return pNtk->GetIntIndex(i) > pNtk->GetIntIndex(j);
+      });
+      break;
+    case 10: // reverse topological order + random PI
+      SetRandPiOrder();
+      pNtk->SortFanins(id, [&](int i, bool ci, int j, bool cj) {
+        if(pNtk->IsPi(i) && pNtk->IsPi(j)) {
+          return vRandPiOrder[pNtk->GetPiIndex(i)] > vRandPiOrder[pNtk->GetPiIndex(j)];
+        }
+        if(pNtk->IsPi(i)) {
+          return false;
+        }
+        if(pNtk->IsPi(j)) {
+          return true;
+        }
+        return pNtk->GetIntIndex(i) > pNtk->GetIntIndex(j);
+      });
+      break;
+    case 11: // topo + fanout + PI
+      pNtk->SortFanins(id, [&](int i, bool ci, int j, bool cj) {
+        if(pNtk->IsPi(i) && !pNtk->IsPi(j)) {
+          return false;
+        }
+        if(!pNtk->IsPi(i) && pNtk->IsPi(j)) {
+          return true;
+        }
+        if(pNtk->GetNumFanouts(i) > pNtk->GetNumFanouts(j)) {
+          return false;
+        }
+        if(pNtk->GetNumFanouts(i) < pNtk->GetNumFanouts(j)) {
+          return true;
+        }
+        if(pNtk->IsPi(i) && pNtk->IsPi(j)) {
+          return false;
+        }
+        return pNtk->GetIntIndex(i) > pNtk->GetIntIndex(j);
+      });
+      break;
+    case 12: // topo + fanout + sorted PI
+      pNtk->SortFanins(id, [&](int i, bool ci, int j, bool cj) {
+        if(pNtk->IsPi(i) && pNtk->IsPi(j)) {
+          return pNtk->GetPiIndex(i) > pNtk->GetPiIndex(j);
+        }
+        if(pNtk->IsPi(i)) {
+          return false;
+        }
+        if(pNtk->IsPi(j)) {
+          return true;
+        }
+        if(pNtk->GetNumFanouts(i) > pNtk->GetNumFanouts(j)) {
+          return false;
+        }
+        if(pNtk->GetNumFanouts(i) < pNtk->GetNumFanouts(j)) {
+          return true;
+        }
+        return pNtk->GetIntIndex(i) > pNtk->GetIntIndex(j);
+      });
+      break;
+    case 13: // topo + fanout + random PI
+      SetRandPiOrder();
+      pNtk->SortFanins(id, [&](int i, bool ci, int j, bool cj) {
+        if(pNtk->IsPi(i) && pNtk->IsPi(j)) {
+          return vRandPiOrder[pNtk->GetPiIndex(i)] > vRandPiOrder[pNtk->GetPiIndex(j)];
+        }
+        if(pNtk->IsPi(i)) {
+          return false;
+        }
+        if(pNtk->IsPi(j)) {
+          return true;
+        }
+        if(pNtk->GetNumFanouts(i) > pNtk->GetNumFanouts(j)) {
+          return false;
+        }
+        if(pNtk->GetNumFanouts(i) < pNtk->GetNumFanouts(j)) {
+          return true;
+        }
+        return pNtk->GetIntIndex(i) > pNtk->GetIntIndex(j);
+      });
+      break;
+    case 14: // random order
+      SetRandCosts();
+      pNtk->SortFanins(id, [&](int i, bool ci, int j, bool cj) {
+        return vRandCosts[i] > vRandCosts[j];
+      });
+      break;
+    case 15: // random + PI
+      SetRandCosts();
+      pNtk->SortFanins(id, [&](int i, bool ci, int j, bool cj) {
+        if(pNtk->IsPi(i) && !pNtk->IsPi(j)) {
+          return false;
+        }
+        if(!pNtk->IsPi(i) && pNtk->IsPi(j)) {
+          return true;
+        }
+        return vRandCosts[i] > vRandCosts[j];
+      });
+      break;
+    case 16: // random + fanout
+      SetRandCosts();
+      pNtk->SortFanins(id, [&](int i, bool ci, int j, bool cj) {
+        if(pNtk->GetNumFanouts(i) > pNtk->GetNumFanouts(j)) {
+          return false;
+        }
+        if(pNtk->GetNumFanouts(i) < pNtk->GetNumFanouts(j)) {
+          return true;
+        }
+        return vRandCosts[i] > vRandCosts[j];
+      });
+      break;
+    case 17: // random + fanout + PI
+      SetRandCosts();
+      pNtk->SortFanins(id, [&](int i, bool ci, int j, bool cj) {
+        if(pNtk->IsPi(i) && !pNtk->IsPi(j)) {
+          return false;
+        }
+        if(!pNtk->IsPi(i) && pNtk->IsPi(j)) {
+          return true;
+        }
+        if(pNtk->GetNumFanouts(i) > pNtk->GetNumFanouts(j)) {
+          return false;
+        }
+        if(pNtk->GetNumFanouts(i) < pNtk->GetNumFanouts(j)) {
+          return true;
+        }
+        return vRandCosts[i] > vRandCosts[j];
+      });
+      break;
+    default:
+      assert(0);
+    }
+  }
+  
+  /* }}} */
+  
+  /* {{{ Reduce fanin */
+  
+  template <typename Ntk, typename Ana>
+  inline bool Optimizer<Ntk, Ana>::ReduceFanin(int id, bool fRemoveUnused) {
+    assert(pNtk->GetNumFanouts(id) > 0);
+    SortFanins(id);
+    bool fRemoved = false;
+    for(int idx = 0; idx < pNtk->GetNumFanins(id); idx++) {
+      // skip fanins that were just added
+      if(mapNewFanins.count(id)) {
+        int fi = pNtk->GetFanin(id, idx);
+        if(mapNewFanins[id].count(fi)) {
+          continue;
+        }
+      }
+      // reduce
+      if(pAna->CheckRedundancy(id, idx)) {
+        int fi = pNtk->GetFanin(id, idx);
+        pNtk->RemoveFanin(id, idx);
+        fRemoved = true;
+        idx--;
+        if(fRemoveUnused && pNtk->GetNumFanouts(fi) == 0) {
+          pNtk->RemoveUnused(fi, true);
+        }
+      }
+    }
+    return fRemoved;
+  }
+
+  template <typename Ntk, typename Ana>
+  inline bool Optimizer<Ntk, Ana>::ReduceFaninOneRandom(int id, bool fRemoveUnused) {
+    assert(pNtk->GetNumFanouts(id) > 0);
+    // generate random order
+    vTmp.resize(pNtk->GetNumFanins(id));
+    std::iota(vTmp.begin(), vTmp.end(), 0);
+    std::shuffle(vTmp.begin(), vTmp.end(), rng);
+    for(int idx: vTmp) {
+      // skip fanins that were just added
+      if(mapNewFanins.count(id)) {
+        int fi = pNtk->GetFanin(id, idx);
+        if(mapNewFanins[id].count(fi)) {
+          continue;
+        }
+      }
+      // reduce
+      if(pAna->CheckRedundancy(id, idx)) {
+        int fi = pNtk->GetFanin(id, idx);
+        pNtk->RemoveFanin(id, idx);
+        if(fRemoveUnused && pNtk->GetNumFanouts(fi) == 0) {
+          pNtk->RemoveUnused(fi, true);
+        }
+        return true;
+      }
+    }
+    return false;
+  }
+  
+  // TODO: add a method to minimize the size of fanins (check singles, pairs, trios, and so on)?
+  
+  /* }}} */
+  
+  /* {{{ Reduce */
+
+  template <typename Ntk, typename Ana>
+  void Optimizer<Ntk, Ana>::Reduce() {
+    std::vector<int> vInts = pNtk->GetInts();
+    for(critr it = vInts.rbegin(); it != vInts.rend(); it++) {
+      if(!pNtk->IsInt(*it)) {
+        continue;
+      }
+      if(pNtk->GetNumFanouts(*it) == 0) {
+        pNtk->RemoveUnused(*it);
+        continue;
+      }
+      ReduceFanin(*it);
+      if(pNtk->GetNumFanins(*it) <= 1) {
+        pNtk->Propagate(*it);
+      }
+      /*
+      if(pNtk->GetNumFanins(*it) == 1) {
+        pNtk->RemoveBuffer(*it);
+      }
+      if(pNtk->GetNumFanins(*it) == 0) {
+        pNtk->RemoveConst(*it);
+      }
+      */
+    }
+  }
+
+  template <typename Ntk, typename Ana>
+  void Optimizer<Ntk, Ana>::ReduceRandom() {
+    pNtk->Sweep(false);
+    std::vector<int> vInts = pNtk->GetInts();
+    std::shuffle(vInts.begin(), vInts.end(), rng);
+    for(citr it = vInts.begin(); it != vInts.end(); it++) {
+      if(!pNtk->IsInt(*it)) {
+        continue;
+      }
+      if(pNtk->GetNumFanouts(*it) == 0) {
+        pNtk->RemoveUnused(*it);
+        continue;
+      }
+      ReduceFanin(*it, true);
+      if(pNtk->GetNumFanins(*it) <= 1) {
+        pNtk->Propagate(*it);
+      }
+    }
+  }
+
+  /* }}} */
+
+  /* {{{ Redundancy removal */
+  
+  template <typename Ntk, typename Ana>
+  void Optimizer<Ntk, Ana>::RemoveRedundancy() {
+    std::vector<int> vInts = pNtk->GetInts();
+    for(critr it = vInts.rbegin(); it != vInts.rend();) {
+      if(!pNtk->IsInt(*it)) {
+        it++;
+        continue;
+      }
+      if(pNtk->GetNumFanouts(*it) == 0) {
+        pNtk->RemoveUnused(*it);
+        it++;
+        continue;
+      }
+      bool fReduced = ReduceFanin(*it);
+      if(pNtk->GetNumFanins(*it) <= 1) {
+        pNtk->Propagate(*it);
+      }
+      if(fReduced) {
+        it = vInts.rbegin();
+      } else {
+        it++;
+      }
+    }
+  }
+
+  template <typename Ntk, typename Ana>
+  void Optimizer<Ntk, Ana>::RemoveRedundancyRandom() {
+    pNtk->Sweep(false);
+    std::vector<int> vInts = pNtk->GetInts();
+    std::shuffle(vInts.begin(), vInts.end(), rng);
+    for(citr it = vInts.begin(); it != vInts.end();) {
+      if(!pNtk->IsInt(*it)) {
+        it++;
+        continue;
+      }
+      if(pNtk->GetNumFanouts(*it) == 0) {
+        pNtk->RemoveUnused(*it);
+        it++;
+        continue;
+      }
+      bool fReduced = ReduceFaninOneRandom(*it, true);
+      if(pNtk->GetNumFanins(*it) <= 1) {
+        pNtk->Propagate(*it);
+      }
+      if(fReduced) {
+        std::shuffle(vInts.begin(), vInts.end(), rng);
+        it = vInts.begin();
+      } else {
+        it++;
+      }
+    }
+  }
+
+  /* }}} */
+
+  /* {{{ Addition */
+
+  template <typename Ntk, typename Ana>
+  template <typename T>
+  T Optimizer<Ntk, Ana>::SingleAdd(int id, T begin, T end) {
+    MarkTfo(id);
+    pNtk->ForEachFanin(id, [&](int fi, bool c) {
+      vMarks[fi] = true;
+    });
+    T it = begin;
+    for(; it != end; it++) {
+      if(!pNtk->IsInt(*it) && !pNtk->IsPi(*it)) {
+        continue;
+      }
+      if(vMarks[*it]) {
+        continue;
+      }
+      if(pAna->CheckFeasibility(id, *it, false)) {
+        pNtk->AddFanin(id, *it, false);
+      } else if(pNtk->UseComplementedEdges() && pAna->CheckFeasibility(id, *it, true)) {
+        pNtk->AddFanin(id, *it, true);
+      } else {
+        continue;
+      }
+      mapNewFanins[id].insert(*it);
+      break;
+    }
+    pNtk->ForEachFanin(id, [&](int fi, bool c) {
+      vMarks[fi] = false;
+    });
+    return it;
+  }
+
+  template <typename Ntk, typename Ana>
+  int Optimizer<Ntk, Ana>::MultiAdd(int id, std::vector<int> const &vCands, int nMax) {
+    MarkTfo(id);
+    pNtk->ForEachFanin(id, [&](int fi, bool c) {
+      vMarks[fi] = true;
+    });
+    int nAdded = 0;
+    for(int cand: vCands) {
+      if(!pNtk->IsInt(cand) && !pNtk->IsPi(cand)) {
+        continue;
+      }
+      if(vMarks[cand]) {
+        continue;
+      }
+      if(pAna->CheckFeasibility(id, cand, false)) {
+        pNtk->AddFanin(id, cand, false);
+      } else if(pNtk->UseComplementedEdges() && pAna->CheckFeasibility(id, cand, true)) {
+        pNtk->AddFanin(id, cand, true);
+      } else {
+        continue;
+      }
+      mapNewFanins[id].insert(cand);
+      nAdded++;
+      if(nAdded == nMax) {
+        break;
+      }
+    }
+    pNtk->ForEachFanin(id, [&](int fi, bool c) {
+      vMarks[fi] = false;
+    });
+    return nAdded;
+  }
+
+  /* }}} */
+  
+  /* {{{ Resub */
+
+  template <typename Ntk, typename Ana>
+  void Optimizer<Ntk, Ana>::SingleResub(int id, std::vector<int> const &vCands, bool fGreedy) {
+    // NOTE: this assumes trivial collapse/decompose does not change cost
+    // let us assume the node is not trivially redundant for now
+    assert(pNtk->GetNumFanouts(id) != 0);
+    assert(pNtk->GetNumFanins(id) > 1);
+    // collapse
+    pNtk->TrivialCollapse(id);
+    // save if wanted
+    int slot;
+    if(fGreedy) {
+      slot = pNtk->Save();
+    }
+    double dCost = CostFunction(pNtk);
+    // main loop
+    for(citr it = vCands.begin(); it != vCands.end(); it++) {
+      if(Timeout()) {
+        break;
+      }
+      if(!pNtk->IsInt(id)) {
+        break;
+      }
+      it = SingleAdd<citr>(id, it, vCands.end());
+      if(it == vCands.end()) {
+        break;
+      }
+      RemoveRedundancy();
+      mapNewFanins.clear();
+      double dNewCost = CostFunction(pNtk);
+      if(nVerbose) {
+        std::cout << "cost: " << dCost << " -> " << dNewCost << std::endl;
+      }
+      if(fGreedy) {
+        if(dNewCost <= dCost) {
+          pNtk->Save(slot);
+          dCost = dNewCost;
+        } else {
+          pNtk->Load(slot);
+        }
+      } else {
+        dCost = dNewCost;
+      }
+    }
+    if(pNtk->IsInt(id)) {
+      pNtk->TrivialDecompose(id);
+    }
+    if(fGreedy) {
+      pNtk->PopBack();
+    }
+  }
+
+  template <typename Ntk, typename Ana>
+  void Optimizer<Ntk, Ana>::MultiResub(int id, std::vector<int> const &vCands, bool fGreedy, int nMax) {
+    // NOTE: this assumes trivial collapse/decompose does not change cost
+    // let us assume the node is not trivially redundant for now
+    assert(pNtk->GetNumFanouts(id) != 0);
+    assert(pNtk->GetNumFanins(id) > 1);
+    // save if wanted
+    int slot;
+    if(fGreedy) {
+      slot = pNtk->Save();
+    }
+    double dCost = CostFunction(pNtk);
+    // collapse
+    pNtk->TrivialCollapse(id);
+    // resub
+    MultiAdd(id, vCands, nMax);
+    RemoveRedundancy();
+    mapNewFanins.clear();
+    RemoveRedundancy();
+    double dNewCost = CostFunction(pNtk);
+    if(nVerbose) {
+      std::cout << "cost: " << dCost << " -> " << dNewCost << std::endl;
+    }
+    if(fGreedy && dNewCost > dCost) {
+      pNtk->Load(slot);
+    }
+    if(pNtk->IsInt(id)) {
+      pNtk->TrivialDecompose(id);
+    }
+    if(fGreedy) {
+      pNtk->PopBack();
+    }
+  }
+  
+  /* }}} */
+
+  /* {{{ Apply */
+
+  template <typename Ntk, typename Ana>
+  void Optimizer<Ntk, Ana>::ApplyReverseTopologically(std::function<void(int)> const &func) {
+    std::vector<int> vInts = pNtk->GetInts();
+    for(critr it = vInts.rbegin(); it != vInts.rend(); it++) {
+      if(Timeout()) {
+        break;
+      }
+      if(!pNtk->IsInt(*it)) {
+        continue;
+      }
+      if(nVerbose) {
+        std::cout << "node " << *it << " (" << std::distance(vInts.crbegin(), it) + 1 << "/" << vInts.size() << ")" << std::endl;
+      }
+      func(*it);
+    }
+  }
+  
+  /* }}} */
+  
+  /* {{{ Constructors */
+  
+  template <typename Ntk, typename Ana>
+  Optimizer<Ntk, Ana>::Optimizer(Ntk *pNtk, Parameter const *pPar, std::function<double(Ntk *)> CostFunction) :
+    pNtk(pNtk),
+    nVerbose(pPar->nOptimizerVerbose),
+    CostFunction(CostFunction),
+    nSortType(pPar->nSortType),
+    nFlow(pPar->nOptimizerFlow),
+    target(-1) {
+    pNtk->AddCallback(std::bind(&Optimizer<Ntk, Ana>::ActionCallback, this, std::placeholders::_1));
+    pAna = new Ana(pNtk, pPar);
+    rng.seed(pPar->iSeed);
+  }
+  
+  template <typename Ntk, typename Ana>
+  Optimizer<Ntk, Ana>::~Optimizer() {
+    delete pAna;
+  }
+  
+  template <typename Ntk, typename Ana>
+  void Optimizer<Ntk, Ana>::UpdateNetwork(Ntk *pNtk_, bool fSame) {
+    pNtk = pNtk_;
+    target = -1;
+    pNtk->AddCallback(std::bind(&Optimizer<Ntk, Ana>::ActionCallback, this, std::placeholders::_1));
+    pAna->UpdateNetwork(pNtk, fSame);
+  }
+  
+  /* }}} */
+
+  /* {{{ Run */
+
+  template <typename Ntk, typename Ana>
+  void Optimizer<Ntk, Ana>::Run(seconds nTimeout_) {
+    nTimeout = nTimeout_;
+    start = GetCurrentTime();
+    switch(nFlow) {
+    case 0:
+      RemoveRedundancy();
+      ApplyReverseTopologically([&](int id) {
+        std::vector<int> vCands = pNtk->GetPisInts();
+        SingleResub(id, vCands);
+      });
+      break;
+    case 1:
+      RemoveRedundancy();
+      ApplyReverseTopologically([&](int id) {
+        std::vector<int> vCands = pNtk->GetInts();
+        MultiResub(id, vCands);
+      });
+      break;
+    case 2: {
+      RemoveRedundancy();
+      double dCost = CostFunction(pNtk);
+      while(true) {
+        ApplyReverseTopologically([&](int id) {
+          std::vector<int> vCands = pNtk->GetPisInts();
+          SingleResub(id, vCands);
+        });
+        ApplyReverseTopologically([&](int id) {
+          std::vector<int> vCands = pNtk->GetInts();
+          MultiResub(id, vCands);
+        });
+        double dNewCost = CostFunction(pNtk);
+        if(dNewCost < dCost) {
+          dCost = dNewCost;
+        } else {
+          break;
+        }
+      }
+      break;
+    }
+    default:
+      assert(0);
+    }
+  }
+
+  template <typename Ntk, typename Ana>
+  void Optimizer<Ntk, Ana>::Randomize() {
+    nSortType = rng() % 18;
+    vRandPiOrder.clear();
+    vRandCosts.clear();
+  }
+  
+  /* }}} */
+  
+}
+
+ABC_NAMESPACE_CXX_HEADER_END
diff --git a/src/opt/rrr/rrrParameter.h b/src/opt/rrr/rrrParameter.h
new file mode 100644
index 000000000..571ded31f
--- /dev/null
+++ b/src/opt/rrr/rrrParameter.h
@@ -0,0 +1,26 @@
+#pragma once
+
+ABC_NAMESPACE_CXX_HEADER_START
+
+namespace rrr {
+
+  struct Parameter {
+    int iSeed = 0;
+    int nWords = 10;
+    int nTimeout = 0;
+    int nSchedulerVerbose = 1;
+    int nOptimizerVerbose = 0;
+    int nAnalyzerVerbose = 0;
+    int nSimulatorVerbose = 0;
+    int nSatSolverVerbose = 0;
+    bool fUseBddCspf = false;
+    bool fUseBddMspf = false;
+    int nConflictLimit = 0;
+    int nSortType = 0;
+    int nOptimizerFlow = 0;
+    int nSchedulerFlow = 0;
+  };
+  
+}
+
+ABC_NAMESPACE_CXX_HEADER_END
diff --git a/src/opt/rrr/rrrSatSolver.h b/src/opt/rrr/rrrSatSolver.h
new file mode 100644
index 000000000..6966eb807
--- /dev/null
+++ b/src/opt/rrr/rrrSatSolver.h
@@ -0,0 +1,546 @@
+#pragma once
+
+#include <iostream>
+#include <iomanip>
+#include <vector>
+
+#include <sat/bsat/satSolver.h>
+
+#include "rrrParameter.h"
+#include "rrrTypes.h"
+
+ABC_NAMESPACE_CXX_HEADER_START
+
+namespace rrr {
+
+  template <typename Ntk>
+  class SatSolver {
+  private:
+    // pointer to network
+    Ntk *pNtk;
+
+    // parameters
+    int nVerbose;
+    int nConflictLimit;
+
+    // data
+    sat_solver *pSat;
+    bool status; // false indicates trivial UNSAT
+    int  target; // node for which miter has been encoded
+    std::vector<int> vVars; // SAT variable for each node
+    std::vector<int> vLits; // temporary storage
+    std::vector<VarValue> vValues; // values in satisfied problem
+    bool fUpdate;
+
+    // statistics
+    int nCalls;
+    int nSats;
+    int nUnsats;
+
+    // callback
+    void ActionCallback(Action const &action);
+
+    // encode
+    void EncodeNode(sat_solver *p, std::vector<int> const &v, int id, int to_negate = -1) const;
+    void EncodeMiter(sat_solver *p, std::vector<int> &v, int id); // create a careset miter where the counterpart has the output of target negated
+    void SetTarget(int id);
+    
+  public:
+    // constructors
+    SatSolver(Ntk *pNtk, Parameter const *pPar);
+    ~SatSolver();
+    void UpdateNetwork(Ntk *pNtk_, bool fSame);
+    
+    // checks
+    SatResult CheckRedundancy(int id, int idx);
+    SatResult CheckFeasibility(int id, int fi, bool c);
+
+    // cex
+    std::vector<VarValue> GetCex();
+  };
+
+  /* {{{ Callback */
+
+  template <typename Ntk>
+  void SatSolver<Ntk>::ActionCallback(Action const &action) {
+    if(target == -1) {
+      return;
+    }
+    switch(action.type) {
+    case REMOVE_FANIN:
+      if(action.id != target) {
+        fUpdate = true;
+      }
+      break;
+    case REMOVE_UNUSED:
+      break;
+    case REMOVE_BUFFER:
+    case REMOVE_CONST:
+      if(action.id == target) {
+        target = -1;
+      }
+      break;
+    case ADD_FANIN:
+      if(action.id != target) {
+        fUpdate = true;
+      }
+      break;
+    case TRIVIAL_COLLAPSE:
+      break;
+    case TRIVIAL_DECOMPOSE:
+      fUpdate = true;
+      break;
+    case SORT_FANINS:
+      break;
+    case SAVE:
+      break;
+    case LOAD:
+      target = -1;
+      break;
+    case POP_BACK:
+      break;
+    default:
+      assert(0);
+    }
+  }
+
+  /* }}} */
+  
+  /* {{{ Encode */
+
+  template <typename Ntk>
+  void SatSolver<Ntk>::EncodeNode(sat_solver *p, std::vector<int> const &v, int id, int to_negate) const {
+    int RetValue;
+    int x = -1, y = -1;
+    bool cx, cy;
+    assert(pNtk->GetNodeType(id) == AND);
+    std::string delim;
+    if(nVerbose) {
+      std::cout << "node " << std::setw(3) << id << ": ";
+    }
+    pNtk->ForEachFanin(id, [&](int fi, bool c) {
+      if(x == -1) {
+        x = v[fi];
+        cx = c ^ (fi == to_negate);
+      } else if(y == -1) {
+        y = v[fi];
+        cy = c ^ (fi == to_negate);
+      } else {
+        int z = sat_solver_addvar(p);
+        if(nVerbose) {
+          std::cout << delim << z << " = " << (cx? "!": "") << x << " & " << (cy? "!": "") << y << std::endl;
+          delim = std::string(10, ' ');
+        }
+        RetValue = sat_solver_add_and(p, z, x, y, cx, cy, 0);
+        assert(RetValue);
+        x = z;
+        cx = false;
+        y = v[fi];
+        cy = c ^ (fi == to_negate);
+      }
+    });
+    if(x == -1) {
+      if(nVerbose) {
+        std::cout << delim << v[id] << " = !0" << std::endl;
+      }
+      RetValue = sat_solver_add_const(p, v[id], 0);
+      assert(RetValue);
+    } else if(y == -1) {
+      if(nVerbose) {
+        std::cout << delim << v[id] << " = " << (cx? "!": "") << x << std::endl;
+      }
+      RetValue = sat_solver_add_buffer(p, v[id], x, cx);
+      assert(RetValue);
+    } else {
+      if(nVerbose) {
+        std::cout << delim << v[id] << " = " << (cx? "!": "") << x << " & " << (cy? "!": "") << y << std::endl;
+      }
+      RetValue = sat_solver_add_and(p, v[id], x, y, cx, cy, 0);
+      assert(RetValue);
+    }
+  }
+  
+  template <typename Ntk>
+  void SatSolver<Ntk>::EncodeMiter(sat_solver *p, std::vector<int> &v, int id) {
+    int RetValue;
+    // reset
+    v.clear();
+    sat_solver_restart(p);
+    status = true;
+    // assign vars for the base
+    int nNodes = pNtk->GetNumNodes();
+    sat_solver_setnvars(p, nNodes);
+    v.resize(nNodes);
+    for(int i = 0; i < nNodes; i++) {
+      v[i] = i;
+    }
+    // constrain const-0
+    RetValue = sat_solver_add_const(p, v[0], 1);
+    assert(RetValue);
+    // encode first circuit
+    if(nVerbose) {
+      std::cout << "encoding network" << std::endl;
+    }
+    pNtk->ForEachInt([&](int id) {
+      EncodeNode(p, v, id);
+    });
+    // always care if it is po
+    if(pNtk->IsPoDriver(id)) {
+      return;
+    }
+    // store po vars (NOTE: it's not a lit)
+    vLits.clear();
+    pNtk->ForEachPoDriver([&](int fi, bool c) {
+      vLits.push_back(v[fi]);
+    });
+    // encode an inverted copy
+    if(nVerbose) {
+      std::cout << "encoding an inverted copy" << std::endl;
+    }
+    pNtk->ForEachTfo(id, false, [&](int fo) {
+      v[fo] = sat_solver_addvar(p);
+      EncodeNode(p, v, fo, id);
+    });
+    // encode miter xors
+    if(nVerbose) {
+      std::cout << "encoding miter xors" << std::endl;
+    }
+    int idx = 0;
+    int n = 0;
+    pNtk->ForEachPoDriver([&](int fi, bool c) {
+      assert(fi != id);
+      if(v[fi] != vLits[idx]) {
+        int x = sat_solver_addvar(p);
+        if(nVerbose) {
+          std::cout << x << " = " << v[fi] << " ^ " << vLits[idx] << std::endl;
+        }
+        RetValue = sat_solver_add_xor(p, x, v[fi], vLits[idx], 0);
+        assert(RetValue);
+        vLits[n] = toLitCond(x, 0);
+        n++;
+      }
+      idx++;
+    });
+    vLits.resize(n);
+    // assign or of xors to 1
+    if(nVerbose) {
+      std::cout << "adding miter output clause" << std::endl;
+      std::cout << "(";
+      std::string delim = "";
+      for(int iLit: vLits) {
+        std::cout << delim << (lit_sign(iLit)? "!": "") << lit_var(iLit);
+        delim = ", ";
+      }
+      std::cout << ")" << std::endl;
+    }
+    if(n == 0) {
+      status = false;
+      return;
+    }
+    RetValue = sat_solver_addclause(p, vLits.data(), vLits.data() + n);
+    assert(RetValue);
+  }
+
+  template <typename Ntk>
+  void SatSolver<Ntk>::SetTarget(int id) {
+    if(!fUpdate && id == target) {
+      return;
+    }
+    fUpdate = false;
+    target = id;
+    EncodeMiter(pSat, vVars, target);
+  }
+
+  /* }}} */
+
+  /* {{{ Constructors */
+
+  template <typename Ntk>
+  SatSolver<Ntk>::SatSolver(Ntk *pNtk, Parameter const *pPar) :
+    pNtk(pNtk),
+    nVerbose(pPar->nSatSolverVerbose),
+    nConflictLimit(pPar->nConflictLimit),
+    pSat(sat_solver_new()),
+    status(false),
+    target(-1),
+    fUpdate(false),
+    nCalls(0),
+    nSats(0),
+    nUnsats(0) {
+    pNtk->AddCallback(std::bind(&SatSolver<Ntk>::ActionCallback, this, std::placeholders::_1));
+  }
+
+  template <typename Ntk>
+  SatSolver<Ntk>::~SatSolver() {
+    sat_solver_delete(pSat);
+    std::cout << "SAT solver stats: calls = " << nCalls << " (SAT = " << nSats << ", UNSAT = " << nUnsats << ", UNDET = " << nCalls - nSats - nUnsats << ")" << std::endl;
+  }
+
+  template <typename Ntk>
+  void SatSolver<Ntk>::UpdateNetwork(Ntk *pNtk_, bool fSame) {
+    pNtk = pNtk_;
+    status = false;
+    target = -1;
+    fUpdate = false;
+    pNtk->AddCallback(std::bind(&SatSolver<Ntk>::ActionCallback, this, std::placeholders::_1));
+  }
+
+  /* }}} */
+
+  /* {{{ Checks */
+  
+  template <typename Ntk>
+  SatResult SatSolver<Ntk>::CheckRedundancy(int id, int idx) {
+    SetTarget(id);
+    if(!status) {
+      if(nVerbose) {
+        std::cout << "trivially UNSATISFIABLE" << std::endl;
+      }
+      return UNSAT;
+    }
+    vLits.clear();
+    assert(pNtk->GetNodeType(id) == AND);    
+    pNtk->ForEachFaninIdx(id, [&](int idx2, int fi, bool c) {
+      if(idx == idx2) {
+        vLits.push_back(toLitCond(vVars[fi], !c));
+      } else {
+        vLits.push_back(toLitCond(vVars[fi], c));
+      }
+    });
+    if(nVerbose) {
+      std::cout << "solving with assumptions: ";
+      std::string delim = "";
+      for(int iLit: vLits) {
+        std::cout << delim << (lit_sign(iLit)? "!": "") << lit_var(iLit);
+        delim = ", ";
+      }
+      std::cout << std::endl;
+    }
+    nCalls++;
+    int res = sat_solver_solve(pSat, vLits.data(), vLits.data() + vLits.size(), nConflictLimit, 0 /*nInsLimit*/, 0 /*nConfLimitGlobal*/, 0 /*nInsLimitGlobal*/);
+    if(res == l_False) {
+      if(nVerbose) {
+        std::cout << "UNSATISFIABLE" << std::endl;
+      }
+      nUnsats++;
+      return UNSAT;
+    }
+    if(res == l_Undef) {
+      if(nVerbose) {
+        std::cout << "UNDETERMINED" << std::endl;
+      }
+      return UNDET;
+    }
+    assert(res == l_True);
+    if(nVerbose) {
+      std::cout << "SATISFIABLE" << std::endl;
+    }
+    nSats++;
+    vValues.clear();
+    vValues.resize(pNtk->GetNumNodes());
+    pNtk->ForEachPi([&](int id) {
+      if(sat_solver_var_value(pSat, vVars[id])) {
+        vValues[id] = TEMP_TRUE;
+      } else {
+        vValues[id] = TEMP_FALSE;
+      }
+    });
+    pNtk->ForEachInt([&](int id) {
+      if(sat_solver_var_value(pSat, vVars[id])) {
+        vValues[id] = TEMP_TRUE;
+      } else {
+        vValues[id] = TEMP_FALSE;
+      }
+    });
+    // required values
+    pNtk->ForEachFaninIdx(id, [&](int idx2, int fi, bool c) {
+      assert((vValues[fi] == TEMP_TRUE) ^ (idx == idx2) ^ c);
+      vValues[fi] = DecideVarValue(vValues[fi]);
+    });
+    return SAT;
+  }
+
+  template <typename Ntk>
+  SatResult SatSolver<Ntk>::CheckFeasibility(int id, int fi, bool c) {
+    SetTarget(id);
+    if(!status) {
+      if(nVerbose) {
+        std::cout << "trivially UNSATISFIABLE" << std::endl;
+      }
+      return UNSAT;
+    }
+    vLits.clear();
+    assert(pNtk->GetNodeType(id) == AND);
+    vLits.push_back(toLit(vVars[id]));
+    vLits.push_back(toLitCond(vVars[fi], !c));
+    if(nVerbose) {
+      std::cout << "solving with assumptions: ";
+      std::string delim = "";
+      for(int iLit: vLits) {
+        std::cout << delim << (lit_sign(iLit)? "!": "") << lit_var(iLit);
+        delim = ", ";
+      }
+      std::cout << std::endl;
+    }
+    nCalls++;
+    int res = sat_solver_solve(pSat, vLits.data(), vLits.data() + vLits.size(), nConflictLimit, 0 /*nInsLimit*/, 0 /*nConfLimitGlobal*/, 0 /*nInsLimitGlobal*/);
+    if(res == l_False) {
+      if(nVerbose) {
+        std::cout << "UNSATISFIABLE" << std::endl;
+      }
+      nUnsats++;
+      return UNSAT;
+    }
+    if(res == l_Undef) {
+      if(nVerbose) {
+        std::cout << "UNDETERMINED" << std::endl;
+      }
+      return UNDET;
+    }
+    assert(res == l_True);
+    if(nVerbose) {
+      std::cout << "SATISFIABLE" << std::endl;
+    }
+    nSats++;
+    vValues.clear();
+    vValues.resize(pNtk->GetNumNodes());
+    pNtk->ForEachPi([&](int id) {
+      if(sat_solver_var_value(pSat, vVars[id])) {
+        vValues[id] = TEMP_TRUE;
+      } else {
+        vValues[id] = TEMP_FALSE;
+      }
+    });
+    pNtk->ForEachInt([&](int id) {
+      if(sat_solver_var_value(pSat, vVars[id])) {
+        vValues[id] = TEMP_TRUE;
+      } else {
+        vValues[id] = TEMP_FALSE;
+      }
+    });
+    // required values
+    assert(vValues[id] == TEMP_TRUE);
+    vValues[id] = DecideVarValue(vValues[id]);
+    assert((vValues[fi] == TEMP_TRUE) ^ !c);
+    vValues[fi] = DecideVarValue(vValues[fi]);
+    return SAT;
+  }
+  
+  /* }}} */
+
+  /* {{{ Cex */
+
+  template <typename Ntk>
+  std::vector<VarValue> SatSolver<Ntk>::GetCex() {
+    if(nVerbose) {
+      std::cout << "cex: ";
+      pNtk->ForEachPi([&](int id) {
+        std::cout << GetVarValueChar(vValues[id]);
+      });
+      std::cout << std::endl;
+    }
+    // reverse simulation
+    pNtk->ForEachIntReverse([&](int id) {
+      switch(pNtk->GetNodeType(id)) {
+      case AND:
+        if(vValues[id] == TRUE) {
+          pNtk->ForEachFanin(id, [&](int fi, bool c) {
+            assert((vValues[fi] == TEMP_TRUE || vValues[fi] == TRUE) ^ c);
+            vValues[fi] = DecideVarValue(vValues[fi]);
+          });
+        } else if(vValues[id] == FALSE) {
+          bool fFound =  false;
+          pNtk->ForEachFanin(id, [&](int fi, bool c) {
+            if(fFound) {
+              return;
+            }
+            if(c) {
+              if(vValues[fi] == TRUE) {
+                fFound = true;
+              }
+            } else {
+              if(vValues[fi] == FALSE) {
+                fFound = true;
+              }
+            }
+          });
+          if(!fFound) {
+            pNtk->ForEachFanin(id, [&](int fi, bool c) {
+              if(fFound) {
+                return;
+              }
+              if(c) {
+                if(vValues[fi] == TEMP_TRUE) {
+                  fFound = true;
+                  vValues[fi] = DecideVarValue(vValues[fi]);
+                }
+              } else {
+                if(vValues[fi] == TEMP_FALSE) {
+                  fFound = true;
+                  vValues[fi] = DecideVarValue(vValues[fi]);
+                }
+              }
+            });
+          }
+        }
+        break;
+      default:
+        assert(0);
+      }
+    });
+    if(nVerbose) {
+      std::cout << "pex: ";
+      pNtk->ForEachPi([&](int id) {
+        std::cout << GetVarValueChar(vValues[id]);
+      });
+      std::cout << std::endl;
+    }
+    // debug
+    pNtk->ForEachInt([&](int id) {
+      switch(pNtk->GetNodeType(id)) {
+      case AND:
+        if(vValues[id] == TRUE) {
+          pNtk->ForEachFanin(id, [&](int fi, bool c) {
+            assert(c || vValues[fi] == TRUE);
+            assert(!c || vValues[fi] == FALSE);
+          });
+        } else if(vValues[id] == FALSE) {
+          bool fFound =  false;
+          pNtk->ForEachFanin(id, [&](int fi, bool c) {
+            if(fFound) {
+              return;
+            }
+            if(c) {
+              if(vValues[fi] == TRUE) {
+                fFound = true;
+              }
+            } else {
+              if(vValues[fi] == FALSE) {
+                fFound = true;
+              }
+            }
+          });
+          assert(fFound);
+        }
+        break;
+      default:
+        assert(0);
+      }
+      });
+    // retrieve partial cex
+    std::vector<VarValue> vPartialCex;
+    pNtk->ForEachPi([&](int id) {
+      if(vValues[id] == TRUE || vValues[id] == FALSE) {
+        vPartialCex.push_back(vValues[id]);
+      } else {
+        vPartialCex.push_back(UNDEF);
+      }
+    });
+    return vPartialCex;
+  }
+
+  /* }}} */
+
+}
+
+ABC_NAMESPACE_CXX_HEADER_END
diff --git a/src/opt/rrr/rrrScheduler.h b/src/opt/rrr/rrrScheduler.h
new file mode 100644
index 000000000..3b7f3457d
--- /dev/null
+++ b/src/opt/rrr/rrrScheduler.h
@@ -0,0 +1,232 @@
+#pragma once
+
+#include <iostream>
+#include <iomanip>
+#include <random>
+
+#include "rrrParameter.h"
+#include "rrrTypes.h"
+
+#include "rrrAbc.h"
+
+ABC_NAMESPACE_CXX_HEADER_START
+
+namespace rrr {
+
+  template <typename Ntk, typename Opt>
+  class Scheduler {
+  private:
+    // pointer to network
+    Ntk *pNtk;
+
+    // parameters
+    int nVerbose;
+    int nFlow;
+
+    // copy of parameter (maybe updated during the run)
+    Parameter Par;
+
+    // data
+    time_point start;
+
+    // time
+    seconds GetRemainingTime();
+
+  public:
+    // constructor
+    Scheduler(Ntk *pNtk, Parameter const *pPar);
+
+    // run
+    void Run();
+  };
+
+  /* {{{ time */
+
+  template <typename Ntk, typename Opt>
+  seconds Scheduler<Ntk, Opt>::GetRemainingTime() {
+    if(Par.nTimeout == 0) {
+      return 0;
+    }
+    time_point current = GetCurrentTime();
+    seconds nTimeout = Par.nTimeout - DurationInSeconds(start, current);
+    if(nTimeout == 0) { // avoid glitch
+      return -1;
+    }
+    return nTimeout;
+  }
+
+  /* }}} */
+
+  /* {{{ Constructor */
+  
+  template <typename Ntk, typename Opt>
+  Scheduler<Ntk, Opt>::Scheduler(Ntk *pNtk, Parameter const *pPar) :
+    pNtk(pNtk),
+    nVerbose(pPar->nSchedulerVerbose),
+    nFlow(pPar->nSchedulerFlow),
+    Par(*pPar) {
+    // TODO: allocations and other preparations should be done here
+  }
+
+  /* }}} */
+
+  /* {{{ Run */
+  
+  // TODO: have multiplie optimizers running on different windows/partitions
+  // TODO: run ABC on those windows or even on the entire network, maybe as a separate class
+  
+  template <typename Ntk, typename Opt>
+  void Scheduler<Ntk, Opt>::Run() {
+    start = GetCurrentTime();
+    // prepare cost function
+    std::function<double(Ntk *)> CostFunction;
+    CostFunction = [](Ntk *pNtk) {
+      int nTwoInputSize = 0;
+      pNtk->ForEachInt([&](int id) {
+        nTwoInputSize += pNtk->GetNumFanins(id) - 1;
+      });
+      return nTwoInputSize;
+    };
+    // start flow
+    Opt opt(pNtk, &Par, CostFunction);
+    switch(nFlow) {
+    case 0:
+      opt.Run(GetRemainingTime());
+      break;
+    case 1: { // transtoch
+      double dCost = CostFunction(pNtk);
+      double dBestCost = dCost;
+      Ntk best(*pNtk);
+      if(nVerbose) {
+        std::cout << "start: cost = " << dCost << std::endl;
+      }
+      for(int i = 0; i < 10; i++) {
+        if(GetRemainingTime() < 0) {
+          break;
+        }
+        if(i != 0) {
+          Abc9Execute(pNtk, "&if -K 6; &mfs; &st");
+          dCost = CostFunction(pNtk);
+          opt.UpdateNetwork(pNtk, true);
+          if(nVerbose) {
+            std::cout << "hop " << std::setw(3) << i << ": cost = " << dCost << std::endl;
+          }
+        }
+        for(int j = 0; true; j++) {
+          if(GetRemainingTime() < 0) {
+            break;
+          }
+          opt.Randomize();
+          opt.Run(GetRemainingTime());
+          Abc9Execute(pNtk, "&dc2");
+          double dNewCost = CostFunction(pNtk);
+          if(nVerbose) {
+            std::cout << "\tite " << std::setw(3) << j << ": cost = " << dNewCost << std::endl;
+          }
+          if(dNewCost < dCost) {
+            dCost = dNewCost;
+            opt.UpdateNetwork(pNtk, true);
+          } else {
+            break;
+          }
+        }
+        if(dCost < dBestCost) {
+          dBestCost = dCost;
+          best = *pNtk;
+          i = 0;
+        }
+      }
+      *pNtk = best;
+      if(nVerbose) {
+        std::cout << "end: cost = " << dBestCost << std::endl;
+      }
+      break;
+    }
+    case 2: { // deep
+      double dCost = CostFunction(pNtk);
+      Ntk best(*pNtk);
+      int n = 0;
+      std::mt19937 rng(Par.iSeed);
+      if(nVerbose) {
+        std::cout << "start: cost = " << dCost << std::endl;
+      }
+      for(int i = 0; i < 1000000; i++) {
+        if(GetRemainingTime() < 0) {
+          break;
+        }
+        // deepsyn
+        int fUseTwo = 0;
+        std::string pCompress2rs = "balance -l; resub -K 6 -l; rewrite -l; resub -K 6 -N 2 -l; refactor -l; resub -K 8 -l; balance -l; resub -K 8 -N 2 -l; rewrite -l; resub -K 10 -l; rewrite -z -l; resub -K 10 -N 2 -l; balance -l; resub -K 12 -l; refactor -z -l; resub -K 12 -N 2 -l; rewrite -z -l; balance -l";
+        unsigned Rand = rng();
+        int fDch = Rand & 1;
+        //int fCom = (Rand >> 1) & 3;
+        int fCom = (Rand >> 1) & 1;
+        int fFx  = (Rand >> 2) & 1;
+        int KLut = fUseTwo ? 2 + (i % 5) : 3 + (i % 4);
+        std::string pComp;
+        if ( fCom == 3 )
+          pComp = "; &put; " + pCompress2rs + "; " + pCompress2rs + "; " + pCompress2rs + "; &get";
+        else if ( fCom == 2 )
+          pComp = "; &put; " + pCompress2rs + "; " + pCompress2rs + "; &get";
+        else if ( fCom == 1 )
+          pComp = "; &put; " + pCompress2rs + "; &get";
+        else if ( fCom == 0 )
+          pComp = "; &dc2";
+        std::string Command = "&dch";
+        if(fDch)
+          Command += " -f";
+        Command += "; &if -a -K " + std::to_string(KLut) + "; &mfs -e -W 20 -L 20";
+        if(fFx)
+          Command += "; &fx; &st";
+        Command += pComp;
+        Abc9Execute(pNtk, Command);
+        if(nVerbose) {
+          std::cout << "ite " << std::setw(6) << i << ": cost = " << CostFunction(pNtk) << std::endl;
+        }
+        // rrr
+        for(int j = 0; j < n; j++) {
+          if(GetRemainingTime() < 0) {
+            break;
+          }
+          opt.Randomize();
+          opt.UpdateNetwork(pNtk, true);
+          opt.Run(GetRemainingTime());
+          if(rng() & 1) {
+            Abc9Execute(pNtk, "&dc2");
+          } else {
+            Abc9Execute(pNtk, "&put; " + pCompress2rs + "; &get");
+          }
+          if(nVerbose) {
+            std::cout << "\trrr " << std::setw(6) << j << ": cost = " << CostFunction(pNtk) << std::endl;
+          }
+        }
+        // eval
+        double dNewCost = CostFunction(pNtk);
+        if(dNewCost < dCost) {
+          dCost = dNewCost;
+          best = *pNtk;
+        } else {
+          n++;
+        }
+      }
+      *pNtk = best;
+      if(nVerbose) {
+        std::cout << "end: cost = " << dCost << std::endl;
+      }
+      break;
+    }
+    default:
+      assert(0);
+    }
+    time_point end = GetCurrentTime();
+    double elapsed_seconds = Duration(start, end);
+    if(nVerbose) {
+      std::cout << "elapsed: " << std::fixed << std::setprecision(3) << elapsed_seconds << "s" << std::endl;
+    }
+  }
+
+  /* }}} */
+
+}
+
+ABC_NAMESPACE_CXX_HEADER_END
diff --git a/src/opt/rrr/rrrSimulator.h b/src/opt/rrr/rrrSimulator.h
new file mode 100644
index 000000000..67a6c7738
--- /dev/null
+++ b/src/opt/rrr/rrrSimulator.h
@@ -0,0 +1,812 @@
+#pragma once
+
+#include <iostream>
+#include <iomanip>
+#include <vector>
+#include <set>
+#include <algorithm>
+#include <random>
+#include <bitset>
+
+#include "rrrParameter.h"
+#include "rrrTypes.h"
+
+ABC_NAMESPACE_CXX_HEADER_START
+
+namespace rrr {
+
+  template <typename Ntk>
+  class Simulator {
+  private:
+    // aliases
+    using word = unsigned long long;
+    using itr = std::vector<word>::iterator;
+    using citr = std::vector<word>::const_iterator;
+    static constexpr word one = 0xffffffffffffffff;
+
+    // pointer to network
+    Ntk *pNtk;
+    
+    // parameters
+    int nVerbose;
+    int nWords;
+
+    // data
+    int target; // node for which the careset has been computed
+    std::vector<word> vValues;
+    std::vector<word> vValues2; // simulation with an inverter
+    std::vector<word> care; // careset
+    std::vector<word> tmp;
+
+    // partial cex
+    int iPivot;
+    std::vector<word> vAssignedStimuli;
+
+    // updates
+    bool fUpdate;
+    std::set<int> sUpdates;
+
+    // backups
+    std::vector<Simulator> vBackups;
+
+    // statistics
+    int nAdds;
+    int nResets;
+    
+    // vector computations
+    void Clear(int n, itr x) const;
+    void Fill(int n, itr x) const;
+    void Copy(int n, itr dst, citr src, bool c) const;
+    void And(int n, itr dst, citr src0, citr src1, bool c0, bool c1) const;
+    void Xor(int n, itr dst, citr src0, citr src1, bool c) const;
+    bool IsZero(int n, citr x) const;
+    bool IsEq(int n, citr x, citr y) const;
+    void Print(int n, citr x) const;
+
+    // callback
+    void ActionCallback(Action const &action);
+
+    // simulation
+    void SimulateNode(std::vector<word> &v, int id, int to_negate = -1);
+    bool ResimulateNode(std::vector<word> &v, int id, int to_negate = -1);
+    void SimulateOneWordNode(std::vector<word> &v, int id, int offset, int to_negate = -1);
+    void Simulate();
+    void Resimulate();
+    void SimulateOneWord(int offset);
+
+    // generate stimuli
+    void GenerateRandomStimuli();
+
+    // careset computation
+    void ComputeCare(int id);
+
+    // save & load
+    void Save(int slot);
+    void Load(int slot);
+
+  public:
+    // constructors
+    Simulator();
+    Simulator(Ntk *pNtk, Parameter const *pPar);
+    ~Simulator();
+    void UpdateNetwork(Ntk *pNtk_, bool fSame);
+
+    // checks
+    bool CheckRedundancy(int id, int idx);
+    bool CheckFeasibility(int id, int fi, bool c);
+
+    // cex
+    void AddCex(std::vector<VarValue> const &vCex);
+  };
+
+
+  /* {{{ Vector computations */
+  
+  template <typename Ntk>
+  inline void Simulator<Ntk>::Clear(int n, itr x) const {
+    std::fill(x, x + n, 0);
+  }
+
+  template <typename Ntk>
+  inline void Simulator<Ntk>::Fill(int n, itr x) const {
+    std::fill(x, x + n, one);
+  }
+  
+  template <typename Ntk>
+  inline void Simulator<Ntk>::Copy(int n, itr dst, citr src, bool c) const {
+    if(!c) {
+      for(int i = 0; i < n; i++, dst++, src++) {
+        *dst = *src;
+      }
+    } else {
+      for(int i = 0; i < n; i++, dst++, src++) {
+        *dst = ~*src;
+      }
+    }
+  }
+  
+  template <typename Ntk>
+  inline void Simulator<Ntk>::And(int n, itr dst, citr src0, citr src1, bool c0, bool c1) const {
+    if(!c0) {
+      if(!c1) {
+        for(int i = 0; i < n; i++, dst++, src0++, src1++) {
+          *dst = *src0 & *src1;
+        }
+      } else {
+        for(int i = 0; i < n; i++, dst++, src0++, src1++) {
+          *dst = *src0 & ~*src1;
+        }
+      }
+    } else {
+      if(!c1) {
+        for(int i = 0; i < n; i++, dst++, src0++, src1++) {
+          *dst = ~*src0 & *src1;
+        }
+      } else {
+        for(int i = 0; i < n; i++, dst++, src0++, src1++) {
+          *dst = ~*src0 & ~*src1;
+        }
+      }
+    }
+  }
+
+  template <typename Ntk>
+  inline void Simulator<Ntk>::Xor(int n, itr dst, citr src0, citr src1, bool c) const {
+    if(!c) {
+      for(int i = 0; i < n; i++, dst++, src0++, src1++) {
+        *dst = *src0 ^ *src1;
+      }
+    } else {
+      for(int i = 0; i < n; i++, dst++, src0++, src1++) {
+        *dst = *src0 ^ ~*src1;
+      }
+    }
+  }
+
+  template <typename Ntk>
+  inline bool Simulator<Ntk>::IsZero(int n, citr x) const {
+    for(int i = 0; i < n; i++, x++) {
+      if(*x) {
+        return false;
+      }
+    }
+    return true;
+  }
+
+  template <typename Ntk>
+  inline bool Simulator<Ntk>::IsEq(int n, citr x, citr y) const {
+    for(int i = 0; i < n; i++, x++, y++) {
+      if(*x != *y) {
+        return false;
+      }
+    }
+    return true;
+  }
+
+  template <typename Ntk>
+  inline void Simulator<Ntk>::Print(int n, citr x) const {
+    std::cout << std::bitset<64>(*x);
+    x++;
+    for(int i = 1; i < n; i++, x++) {
+      std::cout << std::endl << std::string(10, ' ') << std::bitset<64>(*x);
+    }
+  }
+  
+  /* }}} */
+
+  /* {{{ Callback */
+  
+  template <typename Ntk>
+  void Simulator<Ntk>::ActionCallback(Action const &action) {
+    if(target == -1) {
+      return;
+    }
+    switch(action.type) {
+    case REMOVE_FANIN:
+      if(action.id == target) {
+        fUpdate = true;
+      } else {
+        sUpdates.insert(action.id);
+      }
+      break;
+    case REMOVE_UNUSED:
+      break;
+    case REMOVE_BUFFER:
+    case REMOVE_CONST:
+      if(action.id == target) {
+        if(fUpdate) {
+          for(int fo: action.vFanouts) {
+            sUpdates.insert(fo);
+          }
+          fUpdate = false;
+        }
+        target = -1;
+      } else {
+        if(sUpdates.count(action.id)) {
+          sUpdates.erase(action.id);
+          for(int fo: action.vFanouts) {
+            sUpdates.insert(fo);
+          }
+        }
+      }
+      break;
+    case ADD_FANIN:
+      if(action.id == target) {
+        fUpdate = true;
+      } else {
+        sUpdates.insert(action.id);
+      }
+      break;
+    case TRIVIAL_COLLAPSE:
+      break;
+    case TRIVIAL_DECOMPOSE:
+      vValues.resize(nWords * pNtk->GetNumNodes());
+      SimulateNode(vValues, action.fi);
+      break;
+    case SORT_FANINS:
+      break;
+    case SAVE:
+      Save(action.idx);
+      break;
+    case LOAD:
+      Load(action.idx);
+      break;
+    case POP_BACK:
+      // Do nothing: it may be good to keep the word vector allocated
+      break;
+    default:
+      assert(0);
+    }
+  }
+  
+  /* }}} */
+
+  /* {{{ Simulation */
+  
+  template <typename Ntk>
+  void Simulator<Ntk>::SimulateNode(std::vector<word> &v, int id, int to_negate) {
+    itr x = v.end();
+    itr y = v.begin() + id * nWords;
+    bool cx;
+    switch(pNtk->GetNodeType(id)) {
+    case AND:
+      pNtk->ForEachFanin(id, [&](int fi, bool c) {
+        if(x == v.end()) {
+          x = v.begin() + fi * nWords;
+          cx = c ^ (fi == to_negate);
+        } else {
+          And(nWords, y, x, v.begin() + fi * nWords, cx, c ^ (fi == to_negate));
+          x = y;
+          cx = false;
+        }
+      });
+      if(x == v.end()) {
+        Fill(nWords, y);
+      }
+      break;
+    default:
+      assert(0);
+    }
+  }
+      
+  template <typename Ntk>
+  bool Simulator<Ntk>::ResimulateNode(std::vector<word> &v, int id, int to_negate) {
+    itr x = v.end();
+    bool cx;
+    switch(pNtk->GetNodeType(id)) {
+    case AND:
+      pNtk->ForEachFanin(id, [&](int fi, bool c) {
+        if(x == v.end()) {
+          x = v.begin() + fi * nWords;
+          cx = c ^ (fi == to_negate);
+        } else {
+          And(nWords, tmp.begin(), x, v.begin() + fi * nWords, cx, c ^ (fi == to_negate));
+          x = tmp.begin();
+          cx = false;
+        }
+      });
+      if(x == v.end()) {
+        Fill(nWords, tmp.begin());
+      }
+      break;
+    default:
+      assert(0);
+    }
+    itr y = v.begin() + id * nWords;
+    if(IsEq(nWords, y, tmp.begin())) {
+      return false;
+    }
+    Copy(nWords, y, tmp.begin(), false);
+    return true;
+  }
+  
+  template <typename Ntk>
+  void Simulator<Ntk>::SimulateOneWordNode(std::vector<word> &v, int id, int offset, int to_negate) {
+    itr x = v.end();
+    itr y = v.begin() + id * nWords + offset;
+    bool cx;
+    switch(pNtk->GetNodeType(id)) {
+    case AND:
+      pNtk->ForEachFanin(id, [&](int fi, bool c) {
+        if(x == v.end()) {
+          x = v.begin() + fi * nWords + offset;
+          cx = c ^ (fi == to_negate);
+        } else {
+          And(1, y, x, v.begin() + fi * nWords + offset, cx, c ^ (fi == to_negate));
+          x = y;
+          cx = false;
+        }
+      });
+      if(x == v.end()) {
+        Fill(1, y);
+      }
+      break;
+    default:
+      assert(0);
+    }
+  }
+  
+  template <typename Ntk>
+  void Simulator<Ntk>::Simulate() {
+    if(nVerbose) {
+      std::cout << "simulating" << std::endl;
+    }
+    pNtk->ForEachInt([&](int id) {
+      SimulateNode(vValues, id);
+      if(nVerbose) {
+        std::cout << "node " << std::setw(3) << id << ": ";
+        Print(nWords, vValues.begin() + id * nWords);
+        std::cout << std::endl;
+      }
+    });
+  }
+  
+  template <typename Ntk>
+  void Simulator<Ntk>::Resimulate() {
+    if(nVerbose) {
+      std::cout << "resimulating" << std::endl;
+    }
+    pNtk->ForEachTfosUpdate(sUpdates, false, [&](int fo) {
+      bool fUpdated = ResimulateNode(vValues, fo);
+      if(nVerbose) {
+        std::cout << "node " << std::setw(3) << fo << ": ";
+        Print(nWords, vValues.begin() + fo * nWords);
+        std::cout << std::endl;
+      }
+      return fUpdated;
+    });
+    /* alternative version that updates entire TFO
+    pNtk->ForEachTfos(sUpdates, false, [&](int fo) {
+      SimulateNode(vValues, fo);
+      if(nVerbose) {
+        std::cout << "node " << std::setw(3) << fo << ": ";
+        Print(nWords, vValues.begin() + fo * nWords);
+        std::cout << std::endl;
+      }
+    });
+    */
+  }
+
+  template <typename Ntk>
+  void Simulator<Ntk>::SimulateOneWord(int offset) {
+    if(nVerbose) {
+      std::cout << "simulating word " << offset << std::endl;
+    }
+    pNtk->ForEachInt([&](int id) {
+      SimulateOneWordNode(vValues, id, offset);
+      if(nVerbose) {
+        std::cout << "node " << std::setw(3) << id << ": ";
+        Print(1, vValues.begin() + id * nWords + offset);
+        std::cout << std::endl;
+      }
+    });
+  }
+
+  /* }}} */
+
+  /* {{{ Generate stimuli */
+
+  template <typename Ntk>
+  void Simulator<Ntk>::GenerateRandomStimuli() {
+    if(nVerbose) {
+      std::cout << "generating random stimuli" << std::endl;
+    }
+    std::mt19937_64 rng;
+    pNtk->ForEachPi([&](int id) {
+      for(int i = 0; i < nWords; i++) {
+        vValues[id * nWords + i] = rng();
+      }
+      if(nVerbose) {
+        std::cout << "node " << std::setw(3) << id << ": ";
+        Print(nWords, vValues.begin() + id * nWords);
+        std::cout << std::endl;
+      }
+    });
+    Clear(nWords * pNtk->GetNumPis(), vAssignedStimuli.begin());
+  }
+
+  /* }}} */
+
+  /* {{{ Careset computation */
+  
+  template <typename Ntk>
+  void Simulator<Ntk>::ComputeCare(int id) {
+    if(sUpdates.empty() && id == target) {
+      return;
+    }
+    if(fUpdate) {
+      sUpdates.insert(target);
+      fUpdate = false;
+    }
+    if(!sUpdates.empty()) {
+      Resimulate();
+      sUpdates.clear();
+    }
+    target = id;
+    if(nVerbose) {
+      std::cout << "computing careset of " << target << std::endl;
+    }
+    if(pNtk->IsPoDriver(target)) {
+      Fill(nWords, care.begin());
+      if(nVerbose) {
+        std::cout << "care " << std::setw(3) << target << ": ";
+        Print(nWords, care.begin());
+        std::cout << std::endl;
+      }
+      return;
+    }
+    vValues2 = vValues;
+    pNtk->ForEachTfo(target, false, [&](int id) {
+      SimulateNode(vValues2, id, target);
+      if(nVerbose) {
+        std::cout << "node " << std::setw(3) << id << ": ";
+        Print(nWords, vValues2.begin() + id * nWords);
+        std::cout << std::endl;
+      }
+    });
+    /* alternative version that updates only affected TFO
+    pNtk->ForEachTfoUpdate(target, false, [&](int id) {
+      bool fUpdated = ResimulateNode(vValues2, id, target);
+      if(nVerbose) {
+        std::cout << "node " << std::setw(3) << id << ": ";
+        Print(nWords, vValues2.begin() + id * nWords);
+        std::cout << std::endl;
+      }
+      return fUpdated;
+    });
+    */
+    Clear(nWords, care.begin());
+    pNtk->ForEachPoDriver([&](int fi, bool c) {
+      assert(fi != target);
+      for(int i = 0; i < nWords; i++) {
+        care[i] = care[i] | (vValues[fi * nWords + i] ^ vValues2[fi * nWords + i]);
+      }
+    });
+    if(nVerbose) {
+      std::cout << "care " << std::setw(3) << target << ": ";
+      Print(nWords, care.begin());
+      std::cout << std::endl;
+    }
+  }
+  
+  /* }}} */
+
+  /* {{{ Save & load */
+
+  template <typename Ntk>
+  void Simulator<Ntk>::Save(int slot) {
+    if(slot >= vBackups.size()) {
+      vBackups.resize(slot + 1);
+    }
+    vBackups[slot].nWords = nWords;
+    if(sUpdates.empty()) {
+      vBackups[slot].target = target;
+      vBackups[slot].care = care;
+    } else {
+      vBackups[slot].target = -1;
+      vBackups[slot].care = care;
+    }
+    if(fUpdate) {
+      sUpdates.insert(target);
+      fUpdate = false;
+    }
+    if(!sUpdates.empty()) {
+      Resimulate();
+      sUpdates.clear();
+    }
+    vBackups[slot].vValues = vValues;
+    vBackups[slot].iPivot = iPivot;
+    vBackups[slot].vAssignedStimuli = vAssignedStimuli;
+    target = vBackups[slot].target; // assigned to -1 when careset needs updating
+  }
+
+  template <typename Ntk>
+  void Simulator<Ntk>::Load(int slot) {
+    assert(slot < vBackups.size());
+    nWords  = vBackups[slot].nWords;
+    target  = vBackups[slot].target;
+    vValues = vBackups[slot].vValues;
+    care    = vBackups[slot].care;
+    iPivot  = vBackups[slot].iPivot;
+    vAssignedStimuli = vBackups[slot].vAssignedStimuli;
+    fUpdate = false;
+    sUpdates.clear();
+    tmp.resize(nWords);
+  }
+  
+  /* }}} */
+  
+  /* {{{ Constructors */
+  
+  template <typename Ntk>
+  Simulator<Ntk>::Simulator() :
+    pNtk(NULL),
+    nVerbose(0),
+    nWords(0),
+    target(-1),
+    iPivot(0),
+    fUpdate(false),
+    nAdds(0),
+    nResets(0) {
+  }
+  
+  template <typename Ntk>
+  Simulator<Ntk>::Simulator(Ntk *pNtk, Parameter const *pPar) :
+    pNtk(pNtk),
+    nVerbose(pPar->nSimulatorVerbose),
+    nWords(pPar->nWords),
+    target(-1),
+    iPivot(0),
+    fUpdate(false),
+    nAdds(0),
+    nResets(0) {
+    pNtk->AddCallback(std::bind(&Simulator<Ntk>::ActionCallback, this, std::placeholders::_1));
+    vValues.resize(nWords * pNtk->GetNumNodes());
+    care.resize(nWords);
+    tmp.resize(nWords);
+    vAssignedStimuli.resize(nWords * pNtk->GetNumPis());
+    GenerateRandomStimuli();
+    Simulate();
+  }
+
+  template <typename Ntk>
+  Simulator<Ntk>::~Simulator() {
+    if(pNtk) {
+      std::cout << "simulator stats: added CEXs = " << nAdds << ", resets = " << nResets << std::endl;
+    }
+  }
+
+  template <typename Ntk>
+  void Simulator<Ntk>::UpdateNetwork(Ntk *pNtk_, bool fSame) {
+    pNtk = pNtk_;
+    pNtk->AddCallback(std::bind(&Simulator<Ntk>::ActionCallback, this, std::placeholders::_1));
+    vValues.resize(nWords * pNtk->GetNumNodes());
+    target = -1;
+    iPivot = 0;
+    fUpdate = false;
+    sUpdates.clear();
+    if(!fSame) { // reset stimuli if network function changed
+      vAssignedStimuli.clear();
+      vAssignedStimuli.resize(nWords * pNtk->GetNumPis());
+      GenerateRandomStimuli();
+    }
+    Simulate();
+  }
+
+  /* }}} */
+
+  /* {{{ Checks */
+  
+  template <typename Ntk>
+  bool Simulator<Ntk>::CheckRedundancy(int id, int idx) {
+    ComputeCare(id);
+    switch(pNtk->GetNodeType(id)) {
+    case AND: {
+      itr x = vValues.end();
+      bool cx;
+      pNtk->ForEachFaninIdx(id, [&](int idx2, int fi, bool c) {
+        if(idx == idx2) {
+          return;
+        }
+        if(x == vValues.end()) {
+          x = vValues.begin() + fi * nWords;
+          cx = c;
+        } else {
+          And(nWords, tmp.begin(), x, vValues.begin() + fi * nWords, cx, c);
+          x = tmp.begin();
+          cx = false;
+        }
+      });
+      if(x == vValues.end()) {
+        x = care.begin();
+      } else {
+        And(nWords, tmp.begin(), x, care.begin(), cx, false);
+        x = tmp.begin();
+      }
+      int fi = pNtk->GetFanin(id, idx);
+      bool c = pNtk->GetCompl(id, idx);
+      And(nWords, tmp.begin(), x, vValues.begin() + fi * nWords, false, !c);
+      return IsZero(nWords, tmp.begin());
+    }
+    default:
+      assert(0);
+    }
+    return false;
+  }
+
+  template <typename Ntk>
+  bool Simulator<Ntk>::CheckFeasibility(int id, int fi, bool c) {
+    ComputeCare(id);
+    switch(pNtk->GetNodeType(id)) {
+    case AND: {
+      itr x = vValues.end();
+      bool cx;
+      pNtk->ForEachFanin(id, [&](int fi, bool c) {
+        if(x == vValues.end()) {
+          x = vValues.begin() + fi * nWords;
+          cx = c;
+        } else {
+          And(nWords, tmp.begin(), x, vValues.begin() + fi * nWords, cx, c);
+          x = tmp.begin();
+          cx = false;
+        }
+      });
+      if(x == vValues.end()) {
+        x = care.begin();
+      } else {
+        And(nWords, tmp.begin(), x, care.begin(), cx, false);
+        x = tmp.begin();
+      }
+      And(nWords, tmp.begin(), x, vValues.begin() + fi * nWords, false, !c);
+      return IsZero(nWords, tmp.begin());
+    }
+    default:
+      assert(0);
+    }
+    return false;
+  }
+
+  /* }}} */
+
+  /* {{{ Cex */
+
+  template <typename Ntk>
+  void Simulator<Ntk>::AddCex(std::vector<VarValue> const &vCex) {
+    if(nVerbose) {
+      std::cout << "cex: ";
+      for(VarValue c: vCex) {
+        std::cout << GetVarValueChar(c);
+      }
+      std::cout << std::endl;
+    }
+    // record care pi indices
+    assert((int)vCex.size() == pNtk->GetNumPis());
+    std::vector<int> vCarePiIdxs;
+    for(int idx = 0; idx < pNtk->GetNumPis(); idx++) {
+      switch(vCex[idx]) {
+      case TRUE:
+        vCarePiIdxs.push_back(idx);
+        break;
+      case FALSE:
+        vCarePiIdxs.push_back(idx);
+        break;
+      default:
+        break;
+      }
+    }
+    assert(!vCarePiIdxs.empty());
+    // find compatible word
+    int iWord = 0;
+    std::vector<word> vCompatibleBits(1);
+    itr it = vCompatibleBits.begin();
+    for(; iWord < nWords; iWord++) {
+      Fill(1, it);
+      for(int idx: vCarePiIdxs) {
+        int id = pNtk->GetPi(idx);
+        bool c;
+        if(vCex[idx] == TRUE) {
+          c = false;
+        } else {
+          assert(vCex[idx] == FALSE);
+          c = true;
+        }
+        itr x = vValues.begin() + id * nWords + iWord;
+        itr y = vAssignedStimuli.begin() + idx * nWords + iWord;
+        And(1, tmp.begin(), x, y, !c, false);
+        And(1, it, it, tmp.begin(), false, true);
+        if(IsZero(1, it)) {
+          break;
+        }
+      }
+      if(!IsZero(1, it)) {
+        break;
+      }
+    }
+    // find compatible bit
+    int iBit;
+    if(iWord < nWords) {
+      assert(!IsZero(1, it));
+      iBit = 0;
+      while(!((*it >> iBit) & 1)) {
+        iBit++;
+      }
+      if(nVerbose) {
+        std::cout << "fusing into stimulus word " << iWord << " bit " << iBit << std::endl;
+      }
+    } else {
+      // no bits are compatible, so reset at pivot
+      iWord = iPivot / 64;
+      iBit = iPivot % 64;
+      if(nVerbose) {
+        std::cout << "resetting stimulus word " << iWord << " bit " << iBit << std::endl;
+      }
+      word mask = 1ull << iBit;
+      for(int idx = 0; idx < pNtk->GetNumPis(); idx++) {
+        vAssignedStimuli[idx * nWords + iWord] &= ~mask;
+      }
+      iPivot++;
+      if(iPivot == 64 * nWords) {
+        iPivot = 0;
+      }
+      nResets++;
+    }
+    // update stimulus
+    for(int idx: vCarePiIdxs) {
+      int id = pNtk->GetPi(idx);
+      word mask = 1ull << iBit;
+      if(vCex[idx] == TRUE) {
+        vValues[id * nWords + iWord] |= mask;
+      } else {
+        assert(vCex[idx] == FALSE);
+        vValues[id * nWords + iWord] &= ~mask;
+      }
+      vAssignedStimuli[idx * nWords + iWord] |= mask;
+      if(nVerbose) {
+        std::cout << "node " << std::setw(3) << id << ": ";
+        Print(1, vValues.begin() + id * nWords + iWord);
+        std::cout << std::endl;
+        std::cout << "asgn " << std::setw(3) << id << ": ";
+        Print(1, vAssignedStimuli.begin() + idx * nWords + iWord);
+        std::cout << std::endl;
+      }
+    }
+    // simulate
+    SimulateOneWord(iWord);
+    // recompute care with new stimulus
+    if(target != -1 && !pNtk->IsPoDriver(target)) {
+      if(nVerbose) {
+        std::cout << "recomputing careset of " << target << std::endl;
+      }
+      pNtk->ForEachPi([&](int id) {
+        vValues2[id * nWords + iWord] = vValues[id * nWords + iWord];
+      });
+      pNtk->ForEachInt([&](int id) {
+        vValues2[id * nWords + iWord] = vValues[id * nWords + iWord];
+      });
+      pNtk->ForEachTfo(target, false, [&](int id) {
+        SimulateOneWordNode(vValues2, id, iWord, target);
+        if(nVerbose) {
+          std::cout << "node " << std::setw(3) << id << ": ";
+          Print(1, vValues2.begin() + id * nWords + iWord);
+          std::cout << std::endl;
+        }
+      });
+      Clear(1, care.begin() + iWord);
+      pNtk->ForEachPoDriver([&](int fi, bool c) {
+        assert(fi != target);
+        care[iWord] = care[iWord] | (vValues[fi * nWords + iWord] ^ vValues2[fi * nWords + iWord]);
+      });
+      if(nVerbose) {
+        std::cout << "care " << std::setw(3) << target << ": ";
+        Print(1, care.begin() + iWord);
+        std::cout << std::endl;
+      }
+    }
+    nAdds++;
+  }
+  
+  /* }}} */
+  
+}
+
+ABC_NAMESPACE_CXX_HEADER_END
diff --git a/src/opt/rrr/rrrTypes.h b/src/opt/rrr/rrrTypes.h
new file mode 100644
index 000000000..f98c3fd37
--- /dev/null
+++ b/src/opt/rrr/rrrTypes.h
@@ -0,0 +1,207 @@
+#pragma once
+
+#include <iostream>
+#include <vector>
+#include <string>
+#include <chrono>
+#include <cassert>
+
+ABC_NAMESPACE_CXX_HEADER_START
+
+namespace rrr {
+
+  enum NodeType {
+    PI,
+    PO,
+    AND,
+    XOR,
+    LUT
+  };
+
+  enum SatResult {
+    SAT,
+    UNSAT,
+    UNDET
+  };
+
+  enum VarValue: char {
+    UNDEF,
+    TRUE,
+    FALSE,
+    TEMP_TRUE,
+    TEMP_FALSE
+  };
+
+  /* {{{ VarValue functions */
+  
+  static inline VarValue DecideVarValue(VarValue x) {
+    switch(x) {
+    case UNDEF:
+      assert(0);
+    case TRUE:
+      return TRUE;
+    case FALSE:
+      return FALSE;
+    case TEMP_TRUE:
+      return TRUE;
+    case TEMP_FALSE:
+      return FALSE;
+    default:
+      assert(0);
+    }
+  }
+
+  static inline char GetVarValueChar(VarValue x) {
+    switch(x) {
+    case UNDEF:
+      return 'x';
+    case TRUE:
+      return '1';
+    case FALSE:
+      return '0';
+    case TEMP_TRUE:
+      return 't';
+    case TEMP_FALSE:
+      return 'f';
+    default:
+      assert(0);
+    }
+  }
+
+  /* }}} */
+
+  enum ActionType {
+    NONE,
+    REMOVE_FANIN,
+    REMOVE_UNUSED,
+    REMOVE_BUFFER,
+    REMOVE_CONST,
+    ADD_FANIN,
+    TRIVIAL_COLLAPSE,
+    TRIVIAL_DECOMPOSE,
+    SORT_FANINS,
+    SAVE,
+    LOAD,
+    POP_BACK
+  };
+
+  struct Action {
+    ActionType type = NONE;
+    int id = -1;
+    int idx = -1;
+    int fi = -1;
+    bool c = false;
+    std::vector<int> vFanins;
+    std::vector<int> vIndices;
+    std::vector<int> vFanouts;
+  };
+
+  /* {{{ Action functions */
+
+  static inline void PrintAction(Action action) {
+    switch(action.type) {
+    case REMOVE_FANIN:
+      std::cout << "remove fanin";
+      break;
+    case REMOVE_UNUSED:
+      std::cout << "remove unused";
+      break;
+    case REMOVE_BUFFER:
+      std::cout << "remove buffer";
+      break;
+    case REMOVE_CONST:
+      std::cout << "remove const";
+      break;
+    case ADD_FANIN:
+      std::cout << "add fanin";
+      break;
+    case TRIVIAL_COLLAPSE:
+      std::cout << "trivial collapse";
+      break;
+    case TRIVIAL_DECOMPOSE:
+      std::cout << "trivial decompose";
+      break;
+    case SORT_FANINS:
+      std::cout << "sort fanins";
+      break;
+    case SAVE:
+      std::cout << "save";
+      break;
+    case LOAD:
+      std::cout << "load";
+      break;
+    case POP_BACK:
+      std::cout << "pop back";
+      break;
+    default:
+      assert(0);
+    }
+    std::cout << ":";
+    if(action.id != -1) {
+      std::cout << " node " << action.id;
+    }
+    if(action.fi != -1) {
+      std::cout << " fanin " << (action.c? "!": "") << action.fi;
+    }
+    if(action.idx != -1) {
+      std::cout << " index " << action.idx;
+    }
+    std::cout << std::endl;
+    if(!action.vFanins.empty()) {
+      std::cout << "\t" << "fanins: ";
+      std::string delim = "";
+      for(int fi: action.vFanins) {
+        std::cout << delim << fi;
+        delim = ", ";
+      }
+      std::cout << std::endl;
+    }
+    if(!action.vIndices.empty()) {
+      std::cout << "\t" << "indices: ";
+      std::string delim = "";
+      for(int fi: action.vIndices) {
+        std::cout << delim << fi;
+        delim = ", ";
+      }
+      std::cout << std::endl;
+    }
+    if(!action.vFanouts.empty()) {
+      std::cout << "\t" << "fanouts: ";
+      std::string delim = "";
+      for(int fo: action.vFanouts) {
+        std::cout << delim << fo;
+        delim = ", ";
+      }
+      std::cout << std::endl;
+    }
+  }
+
+  /* }}} */
+
+  using seconds = int64_t;
+  using clock_type = std::chrono::steady_clock;
+  using time_point = std::chrono::time_point<clock_type>;
+  
+  /* {{{ Time functions */
+  
+  static inline time_point GetCurrentTime() {
+    return clock_type::now();
+  }
+
+  static inline seconds DurationInSeconds(time_point start, time_point end) {
+    seconds t = (std::chrono::duration_cast<std::chrono::seconds>(end - start)).count();
+    assert(t >= 0);
+    return t;
+  }
+
+  static inline double Duration(time_point start, time_point end) {
+    double t = (std::chrono::duration<double>(end - start)).count();
+    assert(t >= 0);
+    return t;
+  }
+
+  /* }}} */
+  
+}
+
+ABC_NAMESPACE_CXX_HEADER_END