liblts_bisim.h

Go to the documentation of this file.

// Author(s): Muck van Weerdenburg, Jan Friso Groote
// Copyright: see the accompanying file COPYING or copy at
// https://github.com/mCRL2org/mCRL2/blob/master/COPYING
//
// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
//
 
#ifndef _LIBLTS_BISIM_H
#define _LIBLTS_BISIM_H
#include <fstream>
#include "mcrl2/modal_formula/state_formula.h"
#include "mcrl2/lts/lts_utilities.h"
#include "mcrl2/lts/detail/liblts_scc.h"
#include "mcrl2/lts/detail/liblts_merge.h"
#include "mcrl2/lts/lts_aut.h"
#include "mcrl2/lts/lts_fsm.h"
#include "mcrl2/lts/lts_dot.h"
 
namespace mcrl2
{
namespace lts
{
namespace detail
{
 
template < class LTS_TYPE>
class bisim_partitioner
{
 
  public:
    bisim_partitioner(
      LTS_TYPE& l,
      const bool branching=false,
      const bool preserve_divergence=false,
      const bool generate_counter_examples=false)
       : max_state_index(0), 
         aut(l),
         branching(branching),
         store_counter_info(generate_counter_examples)
    {
      assert(branching || !preserve_divergence);
      mCRL2log(log::verbose) << (preserve_divergence?"Divergence preserving b":"B") <<
                  (branching?"ranching b":"") << "isimulation partitioner created for "
                  << l.num_states() << " states and " <<
                  l.num_transitions() << " transitions\n";
 
      if (generate_counter_examples)
      {
        outgoing_transitions = transitions_per_outgoing_state_action_pair(aut.get_transitions(), aut.hidden_label_set());
      }
 
      create_initial_partition(branching,preserve_divergence);
      refine_partition_until_it_becomes_stable(branching, preserve_divergence);
    }
 
 
    ~bisim_partitioner()=default;
 
    void replace_transition_system(const bool branching, const bool preserve_divergences)
    {
      // Put all the non inert transitions in a set. Add the transitions that form a self
      // loop. Such transitions only exist in case divergence preserving branching bisimulation is
      // used. A set is used to remove double occurrences of transitions.
      std::set < transition > resulting_transitions;
 
      for (const transition& i: aut.get_transitions()) 
      {
        const bool is_transition_i_hidden=aut.is_tau(aut.apply_hidden_label_map(i.label()));
        if (!branching ||
            !is_transition_i_hidden ||
            get_eq_class(i.from())!=get_eq_class(i.to()) ||
            (preserve_divergences && i.from()==i.to()))
        {
          resulting_transitions.insert(
            transition(
              get_eq_class(i.from()),
// In the line below all hidden transitions are replaced by an explicit tau. It is possible to 
// always use i.label() where the resulting transition is only implicitly hidden. This has as effect that
// the labels of non inert hidden transitions are preserved (e.g. for use in counter examples) but that
// there are possibly multiple hidden transitions between states, leading to a larger number of 
// transitions than strictly necessary. 
              (is_transition_i_hidden?aut.tau_label_index():i.label()),
              get_eq_class(i.to())));
        }
      }
      // Remove the old transitions
      aut.clear_transitions();
 
      // Copy the transitions from the set into the transition system.
      for (const transition& t: resulting_transitions)
      {
        aut.add_transition(t);
      }
     
      // Merge the states, by setting the state labels of each state to the concatenation of the state labels of its
      // equivalence class. 
      if (aut.has_state_info())   /* If there are no state labels this step can be ignored */
      {
        /* Create a vector for the new labels */
        std::vector<typename LTS_TYPE::state_label_t> new_labels(num_eq_classes());
 
        for(std::size_t i=aut.num_states(); i>0; )
        {
          --i;
          const std::size_t new_index=get_eq_class(i);
          new_labels[new_index] = new_labels[new_index] + aut.state_label(i);
        }
 
        aut.set_num_states(num_eq_classes());
        for(std::size_t i=0; i<num_eq_classes(); ++i)
        {
          aut.set_state_label(i,new_labels[i]);
        }
      }
      else
      {
        aut.set_num_states(num_eq_classes());
      }
      
      aut.set_initial_state(get_eq_class(aut.initial_state()));
    }
 
    std::size_t num_eq_classes() const
    {
      return max_state_index;
    }
 
 
    std::size_t get_eq_class(const std::size_t s) const
    {
      assert(s<block_index_of_a_state.size());
      return blocks[block_index_of_a_state[s]].state_index;
    }
 
 
    bool in_same_class(const std::size_t s, const std::size_t t) const
    {
      return get_eq_class(s)==get_eq_class(t);
    }
 
 
    mcrl2::state_formulas::state_formula counter_formula(const std::size_t s, const std::size_t t);
 
  private:
 
    typedef std::size_t block_index_type;
    typedef std::size_t state_type;
    typedef std::size_t label_type;
 
    state_type max_state_index;
    LTS_TYPE& aut;
    bool branching;
    bool store_counter_info;
 
    struct non_bottom_state
    {
      state_type state;
      std::vector < state_type > inert_transitions; // Only the target state is interesting.
 
      non_bottom_state(const state_type s)
        : state(s)
      {}
      non_bottom_state(const state_type s, const std::vector < state_type > &it)
        : state(s), inert_transitions(it)
      {}
    };
 
    struct block
    {
      state_type state_index;                   // The state number that represent the states in this block
      block_index_type block_index;             // The sequence number of this block.
      block_index_type parent_block_index;      // Index of the parent block.
      // If there is no parent block, this refers to the block
      // itself.
      std::vector < state_type > bottom_states; // The non bottom states must be ordered
      // on tau reachability. The deepest
      // states occur first in the vector.
      std::vector < non_bottom_state > non_bottom_states;
      // The non_inert transitions are grouped per label. The (non-inert) transitions
      // with tau labels are at the end of this vector.
      std::vector < transition > non_inert_transitions;
 
      void swap(block& b)
      {
        state_type state_index1=b.state_index;
        b.state_index=state_index;
        state_index=state_index1;
 
        block_index_type block_index1=b.block_index;
        b.block_index=block_index;
        block_index=block_index1;
 
        block_index_type parent_block_index1=b.parent_block_index;
        b.parent_block_index=parent_block_index;
        parent_block_index=parent_block_index1;
 
        bottom_states.swap(b.bottom_states);
        non_bottom_states.swap(b.non_bottom_states);
        non_inert_transitions.swap(b.non_inert_transitions);
      }
    };
 
    std::vector < block > blocks;
 
    // std::vector < bool > block_is_active;       // Indicates whether this is still a block in the partition.
    // Blocks that are split become inactive.
    std::vector < state_type > block_index_of_a_state;
    std::vector < bool > block_flags;
    std::vector < bool > block_is_in_to_be_processed;
    std::vector < bool > state_flags;
 
    std::vector< block_index_type > to_be_processed;
    std::vector< block_index_type > BL;
 
    // Counter example information based on "Computing Distinguishing Formulas for Branching Bisimulation", 1991 by Henri Korver
    // Given a block B that is split using action alpha and block B' into blocks L and R, with P as the partition just before this split, then
    // right_child[B] = R
    std::map < block_index_type, block_index_type > right_child;
    // split_by_action[B] = alpha,
    std::map < block_index_type, label_type > split_by_action;
    // split_by_block[B] = B',
    std::map < block_index_type, block_index_type > split_by_block;
    // r_alpha[R] = { C \in P | \exists_{s \in R, s' \in C} : s -a-> s'} and
    std::map < block_index_type, std::set < block_index_type > > r_alpha;
    // r_tau[R] = { C \in P | \exists_{s \in R, s' \in C} : s -tau-> s' && C != B}.
    std::map < block_index_type, std::set < block_index_type > > r_tauP;
    // map from source state and action to target state, makes generating counterexample info easier
    outgoing_transitions_per_state_action_t outgoing_transitions;
 
    // A counter for creating fresh variables
    int fresh_var_counter = 0;
 
 
    void create_initial_partition(const bool branching,
                                  const bool preserve_divergences)
 
    {
      to_be_processed.clear();
 
      block initial_partition;
 
      // First store the bottom and non bottom states.
      sort_transitions(aut.get_transitions(), aut.hidden_label_set(), mcrl2::lts::src_lbl_tgt);
 
      state_type last_non_stored_state_number=0;
      bool bottom_state=true;
      std::vector < state_type > current_inert_transitions;
 
      {
        // Reserve enough space, such that no reallocations of the vector are required when adding transitions.
        // For this purpose, first the number of inert transitions must be counted, to avoid reserving too much
        // space. This for instance leads to a waste of memory (terabytes for reducing 30M states), especially,
        // when calculating ia strong bisimulation reduction.
        std::size_t initial_partition_non_inert_counter=0;
        std::size_t current_inert_transition_counter=0;
        const std::vector<transition> & trans=aut.get_transitions();
        for (std::vector<transition>::const_iterator r=trans.begin(); r!=trans.end(); ++r)
        {
          const transition t= *r;
          if (branching && aut.is_tau(aut.apply_hidden_label_map(t.label())))
          {
            if (preserve_divergences && t.from()==t.to())
            {
              initial_partition_non_inert_counter++;
            }
            else
            {
              current_inert_transition_counter++;
            }
          }
        }
        current_inert_transitions.reserve(initial_partition_non_inert_counter);
        initial_partition.non_inert_transitions.reserve(current_inert_transition_counter);
      }
 
      const std::vector<transition>& trans=aut.get_transitions();
      for (std::vector<transition>::const_iterator r=trans.begin(); r!=trans.end(); ++r)
      {
        const transition t=* r;
 
        if (branching && aut.is_tau(aut.apply_hidden_label_map(t.label())))
        {
          if (preserve_divergences && t.from()==t.to())
          {
            initial_partition.non_inert_transitions.push_back(transition(t.from(),aut.tau_label_index(),t.to()));
          }
          else
          {
            current_inert_transitions.push_back(t.to());
            bottom_state=false;
          }
        }
        std::vector<transition>::const_iterator next_i=r;
        ++next_i;
        if (next_i==trans.end() || t.from()!=next_i->from())
        {
          // store the current from state
          for (; last_non_stored_state_number<t.from(); ++last_non_stored_state_number)
          {
            initial_partition.bottom_states.push_back(last_non_stored_state_number);
          }
 
          if (bottom_state)
          {
            initial_partition.bottom_states.push_back(t.from());
          }
          else
          {
            initial_partition.non_bottom_states.push_back(non_bottom_state(t.from()));
            current_inert_transitions.swap(initial_partition.non_bottom_states.back().inert_transitions);
            bottom_state=true;
          }
          assert(last_non_stored_state_number==t.from());
          last_non_stored_state_number++;
        }
      }
 
      for (; last_non_stored_state_number<aut.num_states(); ++last_non_stored_state_number)
      {
        initial_partition.bottom_states.push_back(last_non_stored_state_number);
      }
 
      // Sort the non-bottom states such that the deepest states occur first.
      // Raise an exception if there is a non trivial tau loop.
 
      order_on_tau_reachability(initial_partition.non_bottom_states);
 
      // Store the non-inert transitions (i.e. the non tau transitions)
      sort_transitions(aut.get_transitions(), aut.hidden_label_set(), mcrl2::lts::lbl_tgt_src);
      const std::vector<transition> & trans1=aut.get_transitions();
      for (std::vector<transition>::const_iterator r=trans1.begin(); r!=trans1.end(); ++r)
      {
        const transition t= *r;
        if (!branching || !aut.is_tau(aut.apply_hidden_label_map(t.label())))
        {
          // Note that by sorting the transitions first, the non_inert_transitions are grouped per label.
          initial_partition.non_inert_transitions.push_back(t);
        }
      }
 
      // block_is_active.push_back(true);
      initial_partition.block_index=0;
      initial_partition.state_index=0;
      max_state_index=1;
      initial_partition.parent_block_index=0;
      blocks.push_back(block());
      blocks.back().swap(initial_partition);
      block_index_of_a_state=std::vector < block_index_type >(aut.num_states(),0);
      block_flags.push_back(false);
      state_flags=std::vector < bool >(aut.num_states(),false);
      block_is_in_to_be_processed.push_back(false);
      to_be_processed.clear();
      BL.clear();
    } // end create_initial_partition
 
// Refine the partition until the partition has become stable
#ifndef NDEBUG
    void refine_partition_until_it_becomes_stable(const bool branching, const bool preserve_divergence)
#else
    void refine_partition_until_it_becomes_stable(const bool, const bool)
#endif
    {
#ifndef NDEBUG
      std::size_t consistency_check_counter=1;
      std::size_t consistency_check_barrier=1;
#endif
      bool partition_is_unstable=true; // This boolean indicates that the partition becomes unstable
      // because an inert transition becomes non inert.
      while (!to_be_processed.empty() || partition_is_unstable)
      {
#ifndef NDEBUG
        // Avoid checking too often. This is too time consuming, even in debug mode.
        consistency_check_counter++;
        if (consistency_check_counter>=consistency_check_barrier)
        {
          consistency_check_counter=0;
          consistency_check_barrier=consistency_check_barrier*2;
          check_internal_consistency_of_the_partitioning_data_structure(branching, preserve_divergence);
        }
#endif
        if (to_be_processed.empty() && partition_is_unstable)
        {
          // Put all blocks in to_be_processed;
          for (block_index_type i=0; i< blocks.size(); ++i)
          {
            to_be_processed.push_back(i);
          }
          block_is_in_to_be_processed=std::vector < bool >(blocks.size(),true);
          partition_is_unstable=false;
        }
 
        const block_index_type splitter_index=to_be_processed.back();
        // assert(block_is_in_to_be_processed[splitter_index]||!block_is_active[splitter_index]);
        to_be_processed.pop_back();
        block_is_in_to_be_processed[splitter_index]=false;
 
        // Split with the splitter block, unless it is to_be_processed as we have to reconsider it
        // completely anyhow at some later point.
 
        // Walk through the non_inert_transitions via a natural number, because
        // the blocks and the non_inert_transitions can move around in memory. This makes it
        // unsafe to use interators.
 
        for (std::size_t i=0; i<blocks[splitter_index].non_inert_transitions.size(); ++i)
        {
          // The flag of the starting state of *i is raised and its block is added to BL;
 
          assert(blocks[splitter_index].non_inert_transitions[i].from()<aut.num_states());
          state_flags[blocks[splitter_index].non_inert_transitions[i].from()]=true;
          const block_index_type marked_block_index=block_index_of_a_state[blocks[splitter_index].non_inert_transitions[i].from()];
          if (block_flags[marked_block_index]==false)
          {
            block_flags[marked_block_index]=true;
            BL.push_back(marked_block_index);
          }
 
          // If the label of the next action is different, we must carry out the splitting.
          if ((i+1)==blocks[splitter_index].non_inert_transitions.size() ||
              aut.apply_hidden_label_map(blocks[splitter_index].non_inert_transitions[i].label())!=
              aut.apply_hidden_label_map(blocks[splitter_index].non_inert_transitions[i+1].label()))
          {
            // We consider BL which contains references to all blocks from which a state from splitter
            // can be reached. If not all flags of the non bottom states in a block in BL are set, the
            // non flagged non bottom states are moved to a new block.
 
            split_the_blocks_in_BL(partition_is_unstable,aut.apply_hidden_label_map(blocks[splitter_index].non_inert_transitions[i].label()),splitter_index);
 
          }
 
        }
      }
#ifndef NDEBUG
      check_internal_consistency_of_the_partitioning_data_structure(branching, preserve_divergence);
#endif
      block_flags.clear();
      block_is_in_to_be_processed.clear();
      state_flags.clear();
      to_be_processed.clear();
      BL.clear();
    }
 
 
 
    void refine_partion_with_respect_to_divergences(void);
 
    void split_the_blocks_in_BL(
      bool& partition_is_unstable,
      const label_type splitter_label,
      const block_index_type splitter_block)
    {
      for (std::vector < block_index_type > :: const_iterator i1=BL.begin();
           i1!=BL.end(); ++i1)
      {
        // assert(block_is_active[*i1]);
        block_flags[*i1]=false;
        std::vector < state_type > flagged_states;
        std::vector < state_type > non_flagged_states;
        std::vector < state_type > i1_bottom_states;
        i1_bottom_states.swap(blocks[*i1].bottom_states);
 
        for (std::vector < state_type >::const_iterator j=i1_bottom_states.begin();
             j!=i1_bottom_states.end(); ++j)
        {
          if (state_flags[*j])
          {
            // state is flagged.
            flagged_states.push_back(*j);
          }
          else
          {
            // state is not flagged. It will be moved to a new block.
            non_flagged_states.push_back(*j);
            block_index_of_a_state[*j]=blocks.size();
          }
        }
        assert(!flagged_states.empty()||!blocks[*i1].non_bottom_states.empty()||i1_bottom_states.empty());
        block_index_type reset_state_flags_block=*i1;
 
        if (!non_flagged_states.empty())
        {
          // There are flagged and non flagged states. So, the block must be split.
          // Move the unflagged states to the new block.
 
          if (mCRL2logEnabled(log::debug))
          {
            const std::size_t m = static_cast<std::size_t>(std::pow(10.0, std::floor(std::log10(static_cast<double>((blocks.size()+1)/2)))));
            if ((blocks.size()+1)/2 % m==0)
            {
              mCRL2log(log::debug) << "Bisimulation partitioner: create block " << (blocks.size()+1)/2 << std::endl;
            }
          }
 
          // Create a first new block.
          blocks.push_back(block());
          block_index_type new_block1=blocks.size()-1;
          // block_is_active.push_back(true);
          blocks.back().state_index=max_state_index;
          max_state_index++;
          blocks.back().block_index=new_block1;
          blocks.back().parent_block_index=*i1;
 
          non_flagged_states.swap(blocks.back().bottom_states);
          // Put the indices of first split block to to_be_processed.
          to_be_processed.push_back(blocks.back().block_index);
          block_is_in_to_be_processed.push_back(true);
 
          // Create a second new block.
          blocks.push_back(block());
          block_index_type new_block2=blocks.size()-1;
          // block_is_active.push_back(true);
          blocks.back().state_index=blocks[*i1].state_index;
          blocks.back().block_index=new_block2;
          reset_state_flags_block=new_block2;
          blocks.back().parent_block_index=*i1;
 
          // Move the flagged states to the second block, and let the block index of these states refer to this block.
          flagged_states.swap(blocks.back().bottom_states);
          std::vector < state_type > &reference_to_flagged_states_of_block2=blocks.back().bottom_states;
          for (std::vector < state_type >::const_iterator j=reference_to_flagged_states_of_block2.begin();
               j!=reference_to_flagged_states_of_block2.end(); ++j)
          {
            block_index_of_a_state[*j]=new_block2;
          }
 
          // NOTES: new_block1 = R, new_block2 = L
          // Store counter formula info
          if (store_counter_info)
          {
            right_child[*i1] = new_block1;
            split_by_action[*i1] = splitter_label;
            split_by_block[*i1] = splitter_block;
 
            // compute r_alpha[new_block1]
            std::set<block_index_type> reachable_blocks = {};
            for (state_type source_state : blocks[new_block1].bottom_states)
            {
              for (outgoing_transitions_per_state_action_t::const_iterator i = outgoing_transitions.lower_bound(std::pair<state_type, label_type>(source_state, splitter_label));
                i != outgoing_transitions.upper_bound(std::pair<state_type, label_type>(source_state, splitter_label)); ++i)
              {
                state_type target_state = to(i);
                block_index_type target_block = block_index_of_a_state[target_state];
                reachable_blocks.insert(target_block == new_block1 || target_block == new_block2 ? *i1 : target_block);
              }
            }
            r_alpha[new_block1] = reachable_blocks;
 
            if (branching) {
              // compute r_tau[new_block1]
              std::set<block_index_type> tau_reachable_blocks = {};
              for (state_type source_state : blocks[new_block1].bottom_states)
              {
                for (label_type lab = 0; lab < aut.num_action_labels(); ++lab)
                {
                  if (aut.is_tau(aut.apply_hidden_label_map(lab)))
                  {
                    for (outgoing_transitions_per_state_action_t::const_iterator i = outgoing_transitions.lower_bound(std::pair<state_type, label_type>(source_state, lab));
                      i != outgoing_transitions.upper_bound(std::pair<state_type, label_type>(source_state, lab)); ++i)
                    {
                      state_type target_state = to(i);
                      block_index_type target_block = block_index_of_a_state[target_state];
                      if (!(target_block == *i1 || target_block == new_block1 || target_block == new_block2))
                      {
                        tau_reachable_blocks.insert(target_block);
                      }
                    }
                  }
                }
              }
              r_tauP[new_block1] = tau_reachable_blocks;
            }
          }
 
          // Put the indices of second split block to to_be_processed.
          to_be_processed.push_back(blocks.back().block_index);
          block_is_in_to_be_processed.push_back(true);
 
          // reset the flag of block *i1, which is being split.
          block_is_in_to_be_processed[*i1]=false;
          // block_is_active[*i1]=false;
 
          // The flag fields of the new blocks is set to false;
          block_flags.push_back(false);
          block_flags.push_back(false);
 
          // Declare already some space for transitions, such that we can
          // put inert transitions that become non inert in there when investigating
          // the non bottom states. After investigating the non bottom states,
          // the transitions are split over the vectors below.
 
          std::vector < transition > flagged_non_inert_transitions;
          std::vector < transition > non_flagged_non_inert_transitions;
          // Reserve enough space for transitions to be copied. Otherwise, resizing may lead to
          // lot of unneccesary copying...
          flagged_non_inert_transitions.reserve(blocks[*i1].non_inert_transitions.size());
          non_flagged_non_inert_transitions.reserve(blocks[*i1].non_inert_transitions.size());
 
          // Next we scan the non-bottom states of *i1. If for some non-bottom state the flag is not raised
          // and if none of the outgoing P-inert transitions leads to a state in the old block then this
          // state becomes a non bottom state of B2.
 
          std::vector < non_bottom_state > flagged_non_bottom_states;
          std::vector < non_bottom_state > non_flagged_non_bottom_states;
          std::vector < non_bottom_state > i1_non_bottom_states;
          i1_non_bottom_states.swap(blocks[*i1].non_bottom_states);
          for (typename std::vector < non_bottom_state >::iterator k=i1_non_bottom_states.begin();
               k!=i1_non_bottom_states.end(); ++k)
          {
            const std::vector < state_type > &inert_transitions=k->inert_transitions;
            if (!state_flags[k->state])
            {
              bool all_transitions_end_in_unflagged_block=true;
              for (std::vector < state_type > :: const_iterator l=inert_transitions.begin();
                   all_transitions_end_in_unflagged_block && l!=inert_transitions.end(); ++l)
              {
                if (block_index_of_a_state[*l]!= new_block1)
                {
                  block_index_of_a_state[*l]=new_block2;
                  all_transitions_end_in_unflagged_block=false;
                }
              }
              if (all_transitions_end_in_unflagged_block)
              {
                // Move *k to the non flagged block. Swap the inert transitions to avoid copying.
                non_bottom_state s(k->state);
                s.inert_transitions.swap(k->inert_transitions);
                non_flagged_non_bottom_states.push_back(s);
                block_index_of_a_state[k->state]=new_block1;
                continue;
              }
            }
            // Move *k to the flagged block; note that the transitions can have become
            // non-inert. So, investigate them separately.
            std::vector < state_type > remaining_inert_transitions;
            for (std::vector < state_type > :: const_iterator l=inert_transitions.begin();
                 l!=inert_transitions.end(); ++l)
            {
              if (block_index_of_a_state[*l]==new_block1)
              {
                // The transition *l (*k,tau_label,*l) becomes a non inert transition in the new
                // block.
                non_flagged_non_inert_transitions.push_back(transition(k->state,aut.tau_label_index(),*l));
              }
              else
              {
                // The transition represented by *l remains an inert transition.
                block_index_of_a_state[*l]=new_block2;
                remaining_inert_transitions.push_back(*l);
              }
            }
            if (remaining_inert_transitions.empty()) // The last outgoing inert tau transition just became non inert.
              // k->state has become a bottom state. Otherwise it remains
              // a non bottom state.
            {
              blocks[new_block2].bottom_states.push_back(k->state);
              block_index_of_a_state[k->state]=new_block2;
              partition_is_unstable=true;
            }
            else
            {
              flagged_non_bottom_states.push_back(non_bottom_state(k->state,remaining_inert_transitions));
              block_index_of_a_state[k->state]=new_block2;
            }
          }
          non_flagged_non_bottom_states.swap(blocks[new_block1].non_bottom_states);
          flagged_non_bottom_states.swap(blocks[new_block2].non_bottom_states);
 
          // Finally the non-inert transitions are distributed over both blocks in the obvious way.
          // Note that this must be done after all states are properly put into a new block.
 
          assert(*i1 < blocks.size());
          std::vector < transition > i1_non_inert_transitions;
          i1_non_inert_transitions.swap(blocks[*i1].non_inert_transitions);
          for (std::vector < transition >::iterator k=i1_non_inert_transitions.begin();
               k!=i1_non_inert_transitions.end(); ++k)
          {
            if (block_index_of_a_state[k->to()]==new_block1)
            {
              non_flagged_non_inert_transitions.push_back(*k);
            }
            else
            {
              block_index_of_a_state[k->to()]=new_block2;
              flagged_non_inert_transitions.push_back(*k);
            }
          }
 
          // Only put the parts that we need to store. Avoid that reserved space for
          // current_inert_transition_counter and flagged_non_inert_transitions is stored
          // in the non_inert_transitions in both blocks. Therefore, we do not use a copy
          // constructor, but put elements in there one by one.
          // The two lines below were old code, wasting memory bandwidth.
          // non_flagged_non_inert_transitions.swap(blocks[new_block1].non_inert_transitions);
          // flagged_non_inert_transitions.swap(blocks[new_block2].non_inert_transitions);
 
          blocks[new_block1].non_inert_transitions.reserve(non_flagged_non_inert_transitions.size());
          for(std::vector < transition > ::const_iterator i=non_flagged_non_inert_transitions.begin();
                    i!=non_flagged_non_inert_transitions.end(); i++)
          {
            blocks[new_block1].non_inert_transitions.push_back(*i);
          }
 
          blocks[new_block2].non_inert_transitions.reserve(flagged_non_inert_transitions.size());
          for(std::vector < transition > ::const_iterator i=flagged_non_inert_transitions.begin();
                    i!=flagged_non_inert_transitions.end(); i++)
          {
            blocks[new_block2].non_inert_transitions.push_back(*i);
          }
        }
        else
        {
          // Nothing changed, so put the bottom states back again.
          i1_bottom_states.swap(blocks[*i1].bottom_states);
        }
        // reset the state flags
        std::vector < state_type > &flagged_states_to_be_unflagged=blocks[reset_state_flags_block].bottom_states;
        for (std::vector < state_type >::const_iterator j=flagged_states_to_be_unflagged.begin();
             j!=flagged_states_to_be_unflagged.end(); ++j)
        {
          state_flags[*j]=false;
        }
 
        std::vector < non_bottom_state > &flagged_states_to_be_unflagged1=blocks[reset_state_flags_block].non_bottom_states;
        for (typename std::vector < non_bottom_state >::const_iterator j=flagged_states_to_be_unflagged1.begin();
             j!=flagged_states_to_be_unflagged1.end(); ++j)
        {
          state_flags[j->state]=false;
        }
      }
      BL.clear();
    }
 
    void order_recursively_on_tau_reachability(
      const state_type s,
      std::map < state_type, std::vector < state_type > >& inert_transition_map,
      std::vector < non_bottom_state >& new_non_bottom_states,
      std::set < state_type >& visited)
    {
      if (inert_transition_map.count(s)>0) // The state s is a bottom state. We need not to investigate these.
      {
        if (visited.count(s)==0)
        {
          visited.insert(s);
          std::vector < state_type> &inert_transitions=inert_transition_map[s];
          for (std::vector < state_type>::const_iterator j=inert_transitions.begin();
               j!=inert_transitions.end(); j++)
          {
            order_recursively_on_tau_reachability(*j,inert_transition_map,new_non_bottom_states,visited);
          }
          new_non_bottom_states.push_back(non_bottom_state(s));
          inert_transitions.swap(new_non_bottom_states.back().inert_transitions);
        }
      }
    }
 
    void order_on_tau_reachability(std::vector < non_bottom_state > &non_bottom_states)
    {
      std::set < state_type > visited;
      std::map < state_type, std::vector < state_type > > inert_transition_map;
      for (typename std::vector < non_bottom_state >::iterator i=non_bottom_states.begin();
           i!=non_bottom_states.end(); ++i)
      {
        i->inert_transitions.swap(inert_transition_map[i->state]);
      }
      std::vector < non_bottom_state > new_non_bottom_states;
 
      for (typename std::vector < non_bottom_state >::const_iterator i=non_bottom_states.begin();
           i!=non_bottom_states.end(); ++i)
      {
        order_recursively_on_tau_reachability(i->state, inert_transition_map, new_non_bottom_states, visited);
      }
      new_non_bottom_states.swap(non_bottom_states);
    }
 
    mcrl2::state_formulas::state_formula conjunction(std::set<mcrl2::state_formulas::state_formula> terms) const
    {
      return utilities::detail::join<mcrl2::state_formulas::state_formula>(terms.begin(), terms.end(),
        [](mcrl2::state_formulas::state_formula a, mcrl2::state_formulas::state_formula b)
        { return mcrl2::state_formulas::and_(a, b); }, mcrl2::state_formulas::true_());
    }
 
    regular_formulas::regular_formula create_regular_formula(const mcrl2::lps::multi_action& a) const
    {
      return regular_formulas::regular_formula(action_formulas::multi_action(a.actions()));
    }
 
    regular_formulas::regular_formula create_regular_formula(const mcrl2::lts::action_label_string& a) const
    {
      return mcrl2::regular_formulas::regular_formula(mcrl2::action_formulas::multi_action(
        process::action_list({ process::action(process::action_label(a, {}), {}) })));
    }
 
    mcrl2::state_formulas::state_formula until_formula(const mcrl2::state_formulas::state_formula& phi1, const label_type& a,
                                                       const mcrl2::state_formulas::state_formula& phi2)
    {
      std::string var = "X" + std::to_string(fresh_var_counter++);
      mcrl2::state_formulas::state_formula tauStep =
          mcrl2::state_formulas::may(regular_formulas::regular_formula(action_formulas::multi_action()), mcrl2::state_formulas::variable(var, {}));
      mcrl2::state_formulas::state_formula lastStep = mcrl2::state_formulas::may(create_regular_formula(aut.action_label(a)), phi2);
 
      mcrl2::state_formulas::state_formula until =
          mcrl2::state_formulas::mu(var, {}, mcrl2::state_formulas::and_(phi1, mcrl2::state_formulas::or_(tauStep, lastStep)));
 
      if (aut.is_tau(aut.apply_hidden_label_map(a)))
      {
        until = mcrl2::state_formulas::or_(phi2, until);
      }
      return until;
    }
 
    /* \brief Creates a state formula that distinguishes states in block B1 from state in block B2.
       \details Based on "Computing Distinguishing Formulas for Branching Bisimulation", 1991 by Henri Korver.
                Variable names used below correspond to definitions in this paper.
       \param[in] B1 The block on which the resulting formula must be true.
       \param[in] B2 The block on which the resulting formula must be false.
       \return A distinguishing state formula.
    */
    mcrl2::state_formulas::state_formula counter_formula_aux(const block_index_type B1, const block_index_type B2)
    {
      // First find the smallest block containing both B1 and B2.
      // Find all blocks containing B1.
      std::set < block_index_type > blocks_containing_B1;
      block_index_type B1p = B1;
      blocks_containing_B1.insert(B1p);
      while (blocks[B1p].parent_block_index != B1p)
      {
        B1p = blocks[B1p].parent_block_index;
        blocks_containing_B1.insert(B1p);
      }
 
      // Go up the block tree from B2 until we find a block that contains B1.
      block_index_type DB = B2;
      while (blocks_containing_B1.count(DB) == 0)
      {
        assert(blocks[DB].parent_block_index != DB);
        DB = blocks[DB].parent_block_index;
      }
 
      // Now DB is the smallest block containing both B1 and B2.
 
      const block_index_type R = right_child.at(DB);
      const label_type a = split_by_action.at(DB);
      const block_index_type Bp = split_by_block.at(DB);
 
      // We distinguish Bp with every block in r_alpha[R]
      std::set<mcrl2::state_formulas::state_formula> Gamma2;
      for (block_index_type PPB : r_alpha.at(R))
      {
        Gamma2.insert(counter_formula_aux(Bp, PPB));
      }
      mcrl2::state_formulas::state_formula Phi2 = conjunction(Gamma2);
 
      mcrl2::state_formulas::state_formula Phi;
      if (branching)
      {
        // We distinguish DB with every block in r_tauP[R]
        std::set<mcrl2::state_formulas::state_formula> Gamma1;
        for (block_index_type PPB : r_tauP.at(R))
        {
          Gamma1.insert(counter_formula_aux(DB, PPB));
        }
        mcrl2::state_formulas::state_formula Phi1 = conjunction(Gamma1);
 
        // In case a == tau, we also need to distinguish Bp with R
        if (aut.is_tau(aut.apply_hidden_label_map(a)))
        {
          Phi2 = mcrl2::state_formulas::and_(Phi2, counter_formula_aux(Bp, R));
        }
 
        Phi = until_formula(Phi1, a, Phi2);
      }
      else
      {
        Phi = mcrl2::state_formulas::may(create_regular_formula(aut.action_label(a)), Phi2);
      }
 
      if (blocks_containing_B1.count(R) == 0)
      {
        return Phi;
      }
      else
      {
        return mcrl2::state_formulas::not_(Phi);
      }
    }
 
 
#ifndef NDEBUG
    // The method below is intended to check the consistency of the internal data
    // structure. Its sole purpose is to detect programming errors. It has no
    // side effects on the data structure. If a problem occurs, execution halts with
    // an assert.
 
    void check_internal_consistency_of_the_partitioning_data_structure(
      const bool branching,
      const bool preserve_divergence) const
    {
      state_type total_number_of_states=0;
      std::size_t total_number_of_transitions=0;
 
      assert(!blocks.empty());
      std::set < block_index_type > block_indices;
 
      assert(block_index_of_a_state.size()==aut.num_states());
      for (typename std::vector < block >::const_iterator i=blocks.begin();
           i!=blocks.end(); ++i)
      {
        // Check the block_index.
        assert(i->block_index<blocks.size());
        assert(block_indices.count(i->block_index)==0);
        block_indices.insert(i->block_index);
 
        // Check the bottom states.
        const std::vector < state_type > &i_bottom_states=i->bottom_states;
 
        for (std::vector < state_type >::const_iterator j=i_bottom_states.begin();
             j!=i_bottom_states.end(); ++j)
        {
          total_number_of_states++;
          assert(*j<aut.num_states());
          // Check that the block number of the state is maintained properly.
          assert(block_index_of_a_state[*j]==i->block_index);
        }
 
        // Check the non bottom states. In particular check that there is no tau loop
        // in these non bottom states.
        const std::vector < non_bottom_state > &i_non_bottom_states=i->non_bottom_states;
        std::set < state_type > visited;
        std::set < state_type > local_bottom_states;
        for (typename std::vector < non_bottom_state >::const_iterator j=i_non_bottom_states.begin();
             j!=i_non_bottom_states.end(); ++j)
        {
          local_bottom_states.insert(j->state);
        }
 
        for (typename std::vector < non_bottom_state >::const_iterator j=i_non_bottom_states.begin();
             j!=i_non_bottom_states.end(); ++j)
        {
          total_number_of_states++;
          assert(j->state<aut.num_states());
          // Check that the block number of the state is maintained properly.
          assert(block_index_of_a_state[j->state]==i->block_index);
          const std::vector < state_type > &j_inert_transitions=j->inert_transitions;
          for (std::vector < state_type >::const_iterator k=j_inert_transitions.begin();
               k!=j_inert_transitions.end(); k++)
          {
            total_number_of_transitions++;
            assert(*k<aut.num_states());
            // Check that the inert transitions are well ordered.
            assert(visited.count(*k)>0 || local_bottom_states.count(*k)==0);
          }
          visited.insert(j->state);
        }
 
        // Check the non_inert_transitions. It is required that the transitions
        // are grouped per label, and that tau transitions must be inert.
 
        const std::vector < transition > &i_non_inert_transitions=i->non_inert_transitions;
        std::set < label_type > observed_action_labels;
        for (std::vector < transition >::const_iterator j=i_non_inert_transitions.begin();
             j!=i_non_inert_transitions.end(); ++j)
        {
          total_number_of_transitions++;
          assert(j->to()<aut.num_states());
          assert(j->from()<aut.num_states());
 
          // Check proper grouping of action labels.
          std::vector < transition >::const_iterator j_next=j;
          j_next++;
          if (j_next==i_non_inert_transitions.end() || (aut.apply_hidden_label_map(j->label())!=aut.apply_hidden_label_map(j_next->label())))
          {
            assert(observed_action_labels.count(aut.apply_hidden_label_map(j->label()))==0);
            observed_action_labels.insert(aut.apply_hidden_label_map(j->label()));
          }
 
          // Check whether tau transition in non inert transition vector is inert.
          if (!preserve_divergence && branching && aut.is_tau(aut.apply_hidden_label_map(j->label())))
          {
            assert(j->to()!=j->from());
          }
 
          // Check whether the target state of the transition is in the current block.
          assert(block_index_of_a_state[j->to()]==i->block_index);
        }
      }
 
      // Check total number of states and transitions.
      assert(total_number_of_states==aut.num_states());
      assert(total_number_of_transitions==aut.num_transitions());
 
      // Check block_index_of_a_state
      assert(block_index_of_a_state.size()==aut.num_states());
      for (std::vector < state_type >::const_iterator i=block_index_of_a_state.begin();
           i!=block_index_of_a_state.end(); ++i)
      {
        assert(blocks[*i].block_index== *i);
      }
 
      // Check block_flags that the block flags are all set to false
      for (std::vector < bool >::const_iterator i=block_flags.begin();
           i!=block_flags.end(); ++i)
      {
        assert(!*i);
      }
 
      // Check that state_flags are all false.
      for (std::vector < bool >::const_iterator i=state_flags.begin();
           i!=state_flags.end(); ++i)
      {
        assert(!*i);
      }
 
      // Check to_be_processed
      // Check block_is_in_to_be_processed
      std::vector < bool > temporary_block_is_in_to_be_processed(blocks.size(),false);
 
      for (std::vector< block_index_type > ::const_iterator i=to_be_processed.begin();
           i!=to_be_processed.end(); ++i)
      {
        temporary_block_is_in_to_be_processed[*i]=true;
      }
      for (state_type i=0; i<blocks.size(); ++i)
      {
        assert(!block_is_in_to_be_processed[i] || temporary_block_is_in_to_be_processed[i]);
      }
 
      // Check that BL is empty.
      assert(BL.empty());
 
    }
 
#endif // not NDEBUG
 
};
 
 
template < class LTS_TYPE>
void bisimulation_reduce(
  LTS_TYPE& l,
  const bool branching = false,
  const bool preserve_divergences = false);
 
 
template < class LTS_TYPE>
bool destructive_bisimulation_compare(
  LTS_TYPE& l1,
  LTS_TYPE& l2,
  const bool branching = false,
  const bool preserve_divergences = false,
  const bool generate_counter_examples = false,
  const std::string& counter_example_file = "",
  const bool structured_output = false);
 
 
template < class LTS_TYPE>
bool bisimulation_compare(
  const LTS_TYPE& l1,
  const LTS_TYPE& l2,
  const bool branching = false,
  const bool preserve_divergences = false,
  const bool generate_counter_examples = false,
  const std::string& counter_example_file = "",
  const bool structured_output = false);
 
 
 
 
 
template < class LTS_TYPE>
mcrl2::state_formulas::state_formula bisim_partitioner<LTS_TYPE>::counter_formula(
  const std::size_t s, const std::size_t t)
{
  if (get_eq_class(s)==get_eq_class(t))
  {
    throw mcrl2::runtime_error("Requesting a counter state formula for two bisimilar states. Such a state formula is not useful.");
  }
 
  return counter_formula_aux(block_index_of_a_state[s], block_index_of_a_state[t]);
}
 
 
template < class LTS_TYPE>
void bisimulation_reduce(LTS_TYPE& l,
                         const bool branching /*=false */,
                         const bool preserve_divergences /*=false */)
{
  // First, remove tau loops in case of branching bisimulation.
  if (branching)
  {
    scc_reduce(l,preserve_divergences);
  }
 
  // Secondly, apply the branching bisimulation reduction algorithm. If there are no tau's,
  // this will automatically yield strong bisimulation.
  detail::bisim_partitioner<LTS_TYPE> bisim_part(l, branching, preserve_divergences);
 
  // Assign the reduced LTS
  bisim_part.replace_transition_system(branching,preserve_divergences);
}
 
template < class LTS_TYPE>
bool bisimulation_compare(
  const LTS_TYPE& l1,
  const LTS_TYPE& l2,
  const bool branching /* =false*/,
  const bool preserve_divergences /*=false*/,
  const bool generate_counter_examples /*= false*/,
  const std::string& /* counter_example_file = ""*/,
  const bool structured_output /* = false */)
{
  LTS_TYPE l1_copy(l1);
  LTS_TYPE l2_copy(l2);
  return destructive_bisimulation_compare(l1_copy,l2_copy,branching,preserve_divergences,
                                          generate_counter_examples, structured_output);
}
 
template < class LTS_TYPE>
bool destructive_bisimulation_compare(
  LTS_TYPE& l1,
  LTS_TYPE& l2,
  const bool branching /* =false*/,
  const bool preserve_divergences /*=false*/,
  const bool generate_counter_examples /* = false */,
  const std::string& counter_example_file /*= ""*/,
  const bool /*structured_output = false */)
{
  std::size_t init_l2 = l2.initial_state() + l1.num_states();
  mcrl2::lts::detail::merge(l1,l2);
  l2.clear(); // No use for l2 anymore.
 
  // First remove tau loops in case of branching bisimulation.
  if (branching)
  {
    detail::scc_partitioner<LTS_TYPE> scc_part(l1);
    scc_part.replace_transition_system(preserve_divergences);
    init_l2 = scc_part.get_eq_class(init_l2);
  }
 
  detail::bisim_partitioner<LTS_TYPE> bisim_part(l1, branching, preserve_divergences, generate_counter_examples);
  if (generate_counter_examples && !bisim_part.in_same_class(l1.initial_state(), init_l2))
  {
    mcrl2::state_formulas::state_formula counter_example_formula = bisim_part.counter_formula(l1.initial_state(), init_l2);
 
    std::string filename = "Counterexample.mcf";
    if (!counter_example_file.empty())
    {
      filename = counter_example_file;
    }
    std::ofstream counter_file(filename);
    counter_file << mcrl2::state_formulas::pp(counter_example_formula);
    counter_file.close();
    mCRL2log(mcrl2::log::info) << "Saved counterexample to: \"" << filename << "\"" << std::endl;
 
  }
  return bisim_part.in_same_class(l1.initial_state(),init_l2);
}
 
 
}
}
}
#endif