Merged in rsc-parametric (pull request #1)

PRS -- parametric reaction systems

colour.py

@@ -13,5 +13,16 @@ C_MARK_ERROR = C_BOLD + "[" + C_RED + "!" + C_ENDC + C_BOLD + "]" + C_ENDC
 def colour_str(col, s):
     return col + s + C_ENDC
     
+def green_str(s):
+    return C_GREEN + s + C_ENDC
+    
 def print_error(s):
     print(C_MARK_ERROR + " " + C_RED + s + C_ENDC)
+
+def print_info(s):
+    print("[" + colour_str(C_BOLD, "i") + "] {:s}".format(s))
+
+def print_positive(s):
+    print("[" + colour_str(C_BOLD, "+") + "] {:s}".format(s))
+
+# EOF

logics/bags.py

@@ -1,7 +1,10 @@
 from enum import Enum
-BagDesc_oper = Enum('BagDesc_oper', 'entity true l_and l_or l_not lt le eq ge gt')
+BagDesc_oper = Enum(
+    'BagDesc_oper', 'entity true l_and l_or l_not lt le eq ge gt')
+
 
 class BagDescription(object):
+
     def __init__(self, f_type, L_oper=None, R_oper=None, entity=""):
         self.f_type = f_type
         self.left_operand = L_oper
@@ -20,7 +23,7 @@ class BagDescription(object):
         if self.f_type == BagDesc_oper.l_or:
             return "(" + repr(self.left_operand) + " | " + repr(self.right_operand) + ")"
         if self.f_type == BagDesc_oper.l_not:
-            return "~" + repr(self.left_operand)
+            return "~(" + repr(self.left_operand) + ")"
         if self.f_type == BagDesc_oper.lt:
             return repr(self.left_operand) + " < " + repr(self.right_operand)
         if self.f_type == BagDesc_oper.le:
@@ -34,11 +37,14 @@ class BagDescription(object):
 
     def sanity_check(self):
         """Sanity checks"""
-
         for operand in (self.left_operand, self.right_operand):
             if operand:
-                if not (isinstance(operand, BagDescription) or isinstance(operand, int)):
-                    raise RuntimeError("Unexpected operand type for a bag: " + str(operand))
+                if not(
+                        isinstance(operand, BagDescription)
+                        or isinstance(operand, int)):
+                    raise RuntimeError(
+                        "Unexpected operand type for a bag: {:s} (type: {:s})".format(
+                            str(operand), str(type(operand))))
 
     @classmethod
     def f_entity(cls, entity_name):
@@ -48,6 +54,14 @@ class BagDescription(object):
     def f_TRUE(cls):
         return cls(BagDesc_oper.true)
 
+    @classmethod
+    def f_And(cls, arg_L, arg_R):
+        return cls(BagDesc_oper.l_and, L_oper=arg_L, R_oper=arg_R)
+
+    @classmethod
+    def f_Not(cls, arg):
+        return cls(BagDesc_oper.l_not, L_oper=arg)
+
     def __lt__(self, other):
         return BagDescription(BagDesc_oper.lt, L_oper=self, R_oper=other)
 

logics/rsltl.py

@@ -3,22 +3,34 @@ from enum import Enum
 
 from logics.bags import *
 
-rsLTL_form_type = Enum('rsLTL_form_type', 'bag l_and l_or l_not l_implies t_globally t_finally t_next t_until t_release')
+rsLTL_form_type = Enum(
+    'rsLTL_form_type',
+    'bag l_and l_or l_not l_implies t_globally t_finally t_next t_until t_release')
+
 
 class Formula_rsLTL(object):
 
-    def __init__(self, f_type, L_oper = None, R_oper = None, sub_oper = None, bag = None):
+    def __init__(
+            self, f_type, L_oper=None, R_oper=None, sub_oper=None, bag=None):
         self.f_type = f_type
         self.left_operand = L_oper
         self.right_operand = R_oper
         self.sub_operand = sub_oper
         self.bag_descr = bag
 
+        # it's silly, but it's a frequent mistake
+        if callable(sub_oper):
+            raise RuntimeError(
+                "sub_oper should not be a function. Missing () for {0}?".format(sub_oper))
+
         if isinstance(self.left_operand, BagDescription):
             self.left_operand = Formula_rsLTL.f_bag(self.left_operand)
         if isinstance(self.right_operand, BagDescription):
             self.right_operand = Formula_rsLTL.f_bag(self.right_operand)
             
+        if self.sub_operand is True:
+            self.sub_operand = BagDescription.f_TRUE()
+
     def __repr__(self):
         if self.f_type == rsLTL_form_type.bag:
             return repr(self.bag_descr)
@@ -37,9 +49,15 @@ class Formula_rsLTL(object):
         if self.f_type == rsLTL_form_type.l_implies:
             return "(" + repr(self.left_operand) + " => " + repr(self.right_operand) + ")"
         if self.f_type == rsLTL_form_type.t_until:
-            return "( " + repr(self.left_operand) + " U[" + repr(self.sub_operand) + "] " + repr(self.right_operand) + " )"
+            return "(" + repr(
+                self.left_operand) + " U[" + repr(
+                self.sub_operand) + "] " + repr(
+                self.right_operand) + ")"
         if self.f_type == rsLTL_form_type.t_release:
-            return "( " + repr(self.left_operand) + " R[" + repr(self.sub_operand) + "] " + repr(self.right_operand) + " )"
+            return "(" + repr(
+                self.left_operand) + " R[" + repr(
+                self.sub_operand) + "] " + repr(
+                self.right_operand) + ")"
 
     @property
     def is_bag(self):
@@ -49,6 +67,10 @@ class Formula_rsLTL(object):
     def f_bag(cls, bag_descr):
         return cls(rsLTL_form_type.bag, bag=bag_descr)
 
+    @classmethod
+    def f_Not(cls, arg):
+        return cls(rsLTL_form_type.l_not, L_oper=arg)
+
     @classmethod
     def f_And(cls, arg_L, arg_R):
         return cls(rsLTL_form_type.l_and, L_oper=arg_L, R_oper=arg_R)
@@ -71,7 +93,8 @@ class Formula_rsLTL(object):
 
     @classmethod
     def f_U(cls, sub, arg_L, arg_R):
-        return cls(rsLTL_form_type.t_until, L_oper = arg_L, R_oper = arg_R, sub_oper = sub)
+        return cls(
+            rsLTL_form_type.t_until, L_oper=arg_L, R_oper=arg_R, sub_oper=sub)
 
     @classmethod
     def f_F(cls, sub, arg_L):
@@ -79,7 +102,9 @@ class Formula_rsLTL(object):
 
     @classmethod
     def f_R(cls, sub, arg_L, arg_R):
-        return cls(rsLTL_form_type.t_release, L_oper = arg_L, R_oper = arg_R, sub_oper = sub)
+        return cls(
+            rsLTL_form_type.t_release, L_oper=arg_L, R_oper=arg_R,
+            sub_oper=sub)
 
     def __and__(self, other):
         return Formula_rsLTL(rsLTL_form_type.l_and, L_oper=self, R_oper=other)

logics/rsltl_encoder.py

@@ -8,16 +8,34 @@ class rsLTL_Encoder(object):
         
     def __init__(self, smt_checker):
         self.smt_checker = smt_checker
-        self.v = smt_checker.v
-        self.v_ctx = smt_checker.v_ctx
         self.rs = smt_checker.rs
-        self.loop_position = smt_checker.loop_position
+
+        self.v = None
+        self.v_ctx = None
+        self.loop_position = None
+        
+        # self.load_variables(
+        #     var_rs=smt_checker.v,
+        #     var_ctx=smt_checker.v_ctx,
+        #     var_loop_pos=smt_checker.loop_position)
         
         self.init_ncalls()
         
+    def load_variables(self, var_rs, var_ctx, var_loop_pos):
+        
+        self.v = var_rs
+        self.v_ctx = var_ctx
+        self.loop_position = var_loop_pos
+    
     def get_encoding(self, formula, bound):
+
+        assert self.v is not None
+        assert self.v_ctx is not None
+        assert self.loop_position is not None
+
         self.cache_init(bound)
         self.init_ncalls()
+        
         return self.encode(formula, 0, bound)
         
     def init_ncalls(self):
@@ -103,11 +121,17 @@ class rsLTL_Encoder(object):
         if bag_formula.f_type == BagDesc_oper.gt:
             return self.encode_bag(bag_formula.left_operand, level, context) > int(bag_formula.right_operand)
             
+        assert False, "Unsupported case {:s}".format(bag_formula.f_type)
+    
     def encode_bag_state(self, bag_formula, level):
-        return self.encode_bag(bag_formula, level)
+        res = self.encode_bag(bag_formula, level)
+        assert res is not None
+        return res
         
     def encode_bag_ctx(self, bag_formula, level):
-        return self.encode_bag(bag_formula, level, context=True)
+        res = self.encode_bag(bag_formula, level, context=True)
+        assert res is not None
+        return res
       
     def encode(self, formula, level, bound):
         

rs/__init__.py

@@ -2,6 +2,7 @@ from rs.reaction_system import ReactionSystem
 from rs.context_automaton import ContextAutomaton
 
 from rs.reaction_system_with_concentrations import ReactionSystemWithConcentrations
+from rs.reaction_system_with_concentrations_param import ReactionSystemWithConcentrationsParam, ParameterObj, is_param
 from rs.context_automaton_with_concentrations import ContextAutomatonWithConcentrations
 
 from rs.extended_context_automaton import ExtendedContextAutomaton

rs/reaction_system_with_automaton.py

@@ -1,4 +1,5 @@
 from rs.reaction_system_with_concentrations import ReactionSystemWithConcentrations
+from rs.reaction_system_with_concentrations_param import ReactionSystemWithConcentrationsParam
 from rs.context_automaton_with_concentrations import ContextAutomatonWithConcentrations
 
 class ReactionSystemWithAutomaton(object):
@@ -11,6 +12,17 @@ class ReactionSystemWithAutomaton(object):
         self.rs.show(soft)
         self.ca.show()
         
+    def is_concentr_and_param_compatible(self):
+        """
+        Checks if the underlying RS/CA are compatible 
+        with parameters and concentrations
+        """
+        if not isinstance(self.rs, ReactionSystemWithConcentrationsParam):
+            return False
+        if not isinstance(self.ca, ContextAutomatonWithConcentrations):
+            return False
+        return True
+        
     def is_with_concentrations(self):
         if not isinstance(self.rs, ReactionSystemWithConcentrations):
             return False

rs/reaction_system_with_concentrations.py

@@ -56,8 +56,7 @@ class ReactionSystemWithConcentrations(ReactionSystem):
             if len(e) == 2 and type(e[1]) is int:
                 return True
 
-        print("FATAL. Invalid entity+concentration:")
-        print(e)
+        print("FATAL. Invalid entity+concentration: {:s}".format(e))
         exit(1)
         
         return False

rs/reaction_system_with_concentrations_param.py

@@ -0,0 +1,529 @@
+from sys import exit
+from colour import *
+
+from rs.reaction_system import ReactionSystem
+
+class ParameterObj(object):
+
+    def __init__(self, name):
+        self.name = name
+        
+    def __repr__(self):
+        return "@{0}".format(self.name)
+
+def is_param(some_object):
+    if isinstance(some_object, ParameterObj):
+        return True
+    else:
+        return False
+
+class ReactionSystemWithConcentrationsParam(ReactionSystem):
+
+    def __init__(self):
+
+        self.reactions = []
+        self.parameters = dict()
+        self.meta_reactions = dict()
+        self.permanent_entities = dict()
+        self.background_set = []
+        self.context_entities = []  # legacy. to be removed
+        self.reactions_by_prod = None
+        self.max_concentration = 0
+        self.max_conc_per_ent = dict()
+
+    def add_bg_set_entity(self, e):
+        name = ""
+        def_max_conc = -1
+        if type(e) is tuple and len(e) == 2:
+            name, def_max_conc = e
+        elif type(e) is str:
+            name = e
+            print("\nWARNING: no maximal concentration level specified for:", e, "\n")
+        else:
+            raise RuntimeError(
+                "Bad entity type when adding background set element")
+
+        self.assume_not_in_bgset(name)
+        self.background_set.append(name)
+
+        if def_max_conc != -1:
+            ent_id = self.get_entity_id(name)
+            self.max_conc_per_ent.setdefault(ent_id, 0)
+            if self.max_conc_per_ent[ent_id] < def_max_conc:
+                self.max_conc_per_ent[ent_id] = def_max_conc
+            if self.max_concentration < def_max_conc:
+                self.max_concentration = def_max_conc
+
+    def get_param(self, name):
+        if self.has_param(name):
+            return self.parameters[name]
+        else:
+            param = ParameterObj(name)
+            self.add_param(param)
+            return param
+
+    def has_param(self, name):
+        return name in self.parameters
+
+    def add_param(self, param):
+        if param in self.parameters:
+            raise RuntimeError("Parameter {:s} already exists".format(param))
+        param_key = param.name
+        self.parameters[param_key] = param
+        self.parameters[param_key].idx = len(self.parameters)
+
+    def get_max_concentration_level(self, e):
+
+        if e in self.max_conc_per_ent:
+            return self.max_conc_per_ent[e]
+        else:
+            return self.max_concentration
+
+    def is_valid_entity_with_concentration(self, e):
+        """Sanity check for entities with concentration"""
+
+        if type(e) is tuple:
+            if len(e) == 2 and type(e[1]) is int:
+                return True
+
+        if type(e) is list:
+            if len(e) == 2 and type(e[1]) is int:
+                return True
+
+        print("FATAL. Invalid entity+concentration: {:s}".format(e))
+        exit(1)
+
+        return False
+
+    def get_state_ids(self, state):
+        """Returns entities of the given state without levels"""
+        # return [e for e,c in state]
+        return [self.get_entity_id(e) for e in state]
+
+    def has_non_zero_concentration(self, elem):
+        if elem[1] < 1:
+            raise RuntimeError(
+                "Unexpected concentration level in state: " + str(elem))
+
+    def process_rip(self, R, I, P, ignore_empty_R=False):
+        """Chcecks concentration levels and converts entities names into their ids"""
+
+        if R == [] and not ignore_empty_R:
+            raise RuntimeError("No reactants defined")
+
+        #
+        # REACTANTS
+        #
+        reactants = []
+        if isinstance(R, ParameterObj):
+            reactants = R
+        else:
+            for r in R:
+                self.is_valid_entity_with_concentration(r)
+                self.has_non_zero_concentration(r)
+                entity, level = r
+                reactants.append((self.get_entity_id(entity), level))
+                if self.max_concentration < level:
+                    self.max_concentration = level
+        
+        #
+        # INHIBITORS
+        #
+        inhibitors = []
+        if isinstance(I, ParameterObj):
+            inhibitors = I
+        else:
+            for i in I:
+                self.is_valid_entity_with_concentration(i)
+                self.has_non_zero_concentration(i)
+                entity, level = i
+                inhibitors.append((self.get_entity_id(entity), level))
+                if self.max_concentration < level:
+                    self.max_concentration = level
+        
+        #
+        # PRODUCTS
+        #
+        products = []
+        if isinstance(P, ParameterObj):
+            products = P
+        else:
+            for p in P:
+                self.is_valid_entity_with_concentration(p)
+                self.has_non_zero_concentration(p)
+                entity, level = p
+                products.append((self.get_entity_id(entity), level))
+
+        return reactants, inhibitors, products
+
+    def is_parametric_reaction(self, reaction):
+        result = any([isinstance(r_set, ParameterObj) for r_set in reaction])
+        return result
+
+    def add_reaction(self, R, I, P):
+        """Adds a reaction
+        
+        R, I, and P are sets of entities (not their IDs)
+        """
+        
+        if P == []:
+            raise RuntimeError("No products defined")
+        reaction = self.process_rip(R,I,P)
+        
+        self.reactions.append(reaction)
+
+    def add_reaction_without_reactants(self, R, I, P):
+        """Adds a reaction"""
+        
+        if P == []:
+            raise RuntimeError("No products defined")
+        reaction = self.process_rip(R,I,P,ignore_empty_R=True)
+        self.reactions.append(reaction)
+
+    def add_reaction_inc(self, incr_entity, incrementer, R, I):
+        """Adds a macro/meta reaction for increasing the value of incr_entity"""
+
+        reactants,inhibitors,products = self.process_rip(R,I,[],ignore_empty_R=True)
+        incr_entity_id = self.get_entity_id(incr_entity)
+        self.meta_reactions.setdefault(incr_entity_id,[])
+        self.meta_reactions[incr_entity_id].append(("inc", self.get_entity_id(incrementer), reactants, inhibitors))
+
+    def add_reaction_dec(self, decr_entity, decrementer, R, I):
+        """Adds a macro/meta reaction for decreasing the value of incr_entity"""
+
+        reactants,inhibitors,products = self.process_rip(R,I,[],ignore_empty_R=True)
+        decr_entity_id = self.get_entity_id(decr_entity)
+        self.meta_reactions.setdefault(decr_entity_id,[])
+        self.meta_reactions[decr_entity_id].append(("dec", self.get_entity_id(decrementer), reactants, inhibitors))
+
+    def add_permanency(self, ent, I):
+        """Sets entity to be permanent unless it is inhibited"""
+    
+        ent_id = self.get_entity_id(ent)
+        
+        if ent_id in self.permanent_entities:
+            raise RuntimeError("Permanency for {0} already defined.".format(ent))
+        
+        inhibitors = self.process_rip([],I,[],ignore_empty_R=True)[1]
+        self.permanent_entities[ent_id] = inhibitors
+    
+    def set_context_entities(self, entities):
+        raise NotImplementedError
+
+    def entities_names_set_to_str(self, entities):
+        s = ""
+        for entity in entities:
+            s += entity + ", "
+        s = s[:-2]
+        return s
+
+    def entities_ids_set_to_str(self, entities):
+        s = ""
+        for entity in entities:
+            s += self.get_entity_name(entity) + ", "
+        s = s[:-2]
+        return s
+
+    def state_to_str(self, state):
+        """
+        If state is a parameter, we return
+        the string representation of the whole state
+        which should be the name of the parameter
+        """
+        if isinstance(state, ParameterObj):
+            return str(state)
+        else:
+            s = ""
+            for ent,level in state:
+                s += self.get_entity_name(ent) + "=" + str(level) + ", "
+            s = s[:-2]
+        return s        
+
+    def show_background_set(self):
+        print(C_MARK_INFO + " Background set: {" + self.entities_names_set_to_str(self.background_set) + "}")
+
+    def show_meta_reactions(self):
+        print(C_MARK_INFO + " Meta reactions:")
+        for param_ent,reactions in self.meta_reactions.items():
+            for r_type,command,reactants,inhibitors in reactions:
+                if r_type == "inc" or r_type == "dec":
+                    print(" - [ Type=" + repr(r_type) + " Operand=( " + self.get_entity_name(param_ent) + \
+                        " ) Command=( " + self.get_entity_name(command) + " ) ] -- ( R={" + self.state_to_str(reactants) + "}, I={" + self.state_to_str(inhibitors) + "} )")
+                else:
+                    raise RuntimeError("Unknown meta-reaction type: " + repr(r_type))
+
+    def show_max_concentrations(self):
+        print(C_MARK_INFO + " Maximal allowed concentration levels:")
+        for e,max_conc in self.max_conc_per_ent.items():
+            print(" - {0:^20} = {1:<6}".format(self.get_entity_name(e),max_conc))
+
+    def show_permanent_entities(self):
+        print(C_MARK_INFO + " Permanent entities:")
+        for e,inhibitors in self.permanent_entities.items():
+            print(" - {0:^20}{1:<6}".format(self.get_entity_name(e) + ": ","I={" + self.state_to_str(inhibitors) + "}"))
+
+    def show(self, soft=False):
+        self.show_background_set()
+        self.show_reactions(soft)
+        # self.show_param_reactions(soft)
+        self.show_permanent_entities()
+        self.show_meta_reactions()
+        self.show_max_concentrations()
+        
+    def get_producible_entities(self):
+        """
+        Returns the set of entities that appear as products of 
+        reactions.
+        """
+        
+        producible_entities = set()
+
+        for reaction in self.reactions:
+            product_entities = [e for e,c in reaction[2] if c > 0]
+            producible_entities = producible_entities.union(set(product_entities))
+            
+        return producible_entities
+        
+    def get_reactions_by_product(self):
+        """Sorts reactions by their products and returns a dictionary of products"""
+
+        # assert False
+
+        if self.reactions_by_prod != None:
+            return self.reactions_by_prod
+
+        producible_entities = self.get_producible_entities()
+
+        reactions_by_prod = {}
+
+        for p_e in producible_entities:
+            reactions_by_prod[p_e] = []
+            rcts_for_p_e = reactions_by_prod[p_e]
+
+            for r in self.reactions:
+                product_entities = [e for e,c in r[2]]
+ 
+                if p_e in product_entities:
+                    reactants = r[0]
+                    inhibitors = r[1]
+                    products = [(e,c) for e,c in r[2] if e == p_e]
+                    
+                    prod_conc = products[0][1]
+                    insert_place = None
+                    
+                    # we need to order the reactions w.r.t. the concentration levels produced (increasing order)
+                    for i in range(0,len(rcts_for_p_e)):
+                        
+                        checked_conc = rcts_for_p_e[i][2][0][1]
+                        if prod_conc <= checked_conc:
+                            insert_place = i
+                            break
+
+                    if insert_place == None: # empty or the is only one element which is smaller than the element being added
+                        rcts_for_p_e.append((reactants, inhibitors, products)) # we append (to the end)
+                    else:
+                        rcts_for_p_e.insert(insert_place,(reactants, inhibitors, products))
+                    
+
+        # save in cache
+        self.reactions_by_prod = reactions_by_prod
+
+        return reactions_by_prod
+
+    def get_reaction_system(self):
+        """
+        Translates RSC into RS
+        """
+
+        rs = ReactionSystem()
+
+        for reactants, inhibitors, products in self.reactions:
+
+            new_reactants = []
+            new_inhibitors = []
+            new_products = []
+
+            for ent, conc in reactants:
+                n = self.get_entity_name(ent) + "#" + str(conc)
+                rs.ensure_bg_set_entity(n)
+                new_reactants.append(n)
+
+            for ent, conc in inhibitors:
+                n = self.get_entity_name(ent) + "#" + str(conc)
+                rs.ensure_bg_set_entity(n)
+                new_inhibitors.append(n)
+
+            for ent, conc in products:
+                for i in range(1, conc + 1):
+                    n = self.get_entity_name(ent) + "#" + str(i)
+                    rs.ensure_bg_set_entity(n)
+                    new_products.append(n)
+
+            rs.add_reaction(new_reactants, new_inhibitors, new_products)
+
+        for param_ent, reactions in self.meta_reactions.items():
+            for r_type, command, reactants, inhibitors in reactions:
+
+                param_ent_name = self.get_entity_name(param_ent)
+
+                new_reactants = []
+                new_inhibitors = []
+
+                for ent, conc in reactants:
+                    n = self.get_entity_name(ent) + "#" + str(conc)
+                    rs.ensure_bg_set_entity(n)
+                    new_reactants.append(n)
+
+                for ent, conc in inhibitors:
+                    n = self.get_entity_name(ent) + "#" + str(conc)
+                    rs.ensure_bg_set_entity(n)
+                    new_inhibitors.append(n)
+
+                max_cmd_c = self.max_concentration
+                if command in self.max_conc_per_ent:
+                    max_cmd_c = self.max_conc_per_ent[command]
+                else:
+                    print(
+                        "WARNING:\n\tThere is no maximal concentration level defined for "
+                        + self.get_entity_name(command))
+                    print("\tThis is a very bad idea -- expect degraded performance\n")
+
+                for l in range(1, max_cmd_c + 1):
+
+                    cmd_ent = self.get_entity_name(command) + "#" + str(l)
+                    rs.ensure_bg_set_entity(cmd_ent)
+
+                    if r_type == "inc":
+
+                        # pre_conc  -- predecessor concentration
+                        # succ_conc -- successor concentration concentration
+
+                        for i in range(1, self.max_concentration):
+                            pre_conc = param_ent_name + "#" + str(i)
+                            rs.ensure_bg_set_entity(pre_conc)
+                            new_products = []
+                            succ_value = i + l
+                            for j in range(1, succ_value + 1):
+                                if j > self.max_concentration:
+                                    break
+                                new_p = param_ent_name + "#" + str(j)
+                                rs.ensure_bg_set_entity(new_p)
+                                new_products.append(new_p)
+                            if new_products != []:
+                                rs.add_reaction(
+                                    set(new_reactants + [pre_conc, cmd_ent]),
+                                    set(new_inhibitors),
+                                    set(new_products))
+
+                    elif r_type == "dec":
+                        for i in range(1, self.max_concentration + 1):
+                            pre_conc = param_ent_name + "#" + str(i)
+                            rs.ensure_bg_set_entity(pre_conc)
+                            new_products = []
+                            succ_value = i - l
+                            for j in range(1, succ_value + 1):
+                                if j > self.max_concentration:
+                                    break
+                                new_p = param_ent_name + "#" + str(j)
+                                rs.ensure_bg_set_entity(new_p)
+                                new_products.append(new_p)
+                            if new_products != []:
+                                rs.add_reaction(
+                                    set(new_reactants + [pre_conc, cmd_ent]),
+                                    set(new_inhibitors),
+                                    set(new_products))
+
+                    else:
+                        raise RuntimeError(
+                            "Unknown meta-reaction type: " + repr(r_type))
+
+        for ent, inhibitors in self.permanent_entities.items():
+
+            max_c = self.max_concentration
+            if ent in self.max_conc_per_ent:
+                max_c = self.max_conc_per_ent[ent]
+            else:
+                print(
+                    "WARNING:\n\tThere is no maximal concentration level defined for "
+                    + self.get_entity_name(ent))
+                print("\tThis is a very bad idea -- expect degraded performance\n")
+
+            def e_value(i):
+                return self.get_entity_name(ent) + "#" + str(i)
+
+            for value in range(1, max_c + 1):
+
+                new_reactants = []
+                new_inhibitors = []
+                new_products = []
+
+                new_reactants = [e_value(value)]
+
+                for e_inh, conc in inhibitors:
+                    n = self.get_entity_name(e_inh) + "#" + str(conc)
+                    rs.ensure_bg_set_entity(n)
+                    new_inhibitors.append(n)
+
+                for i in range(1, value + 1):
+                    new_products.append(e_value(i))
+
+                rs.add_reaction(new_reactants, new_inhibitors, new_products)
+
+        return rs
+
+
+# class ReactionSystemWithAutomaton(object):
+#
+#     def __init__(self, reaction_system, context_automaton):
+#         self.rs = reaction_system
+#         self.ca = context_automaton
+#
+#     def show(self, soft=False):
+#         self.rs.show(soft)
+#         self.ca.show()
+#
+#     def is_concentr_and_param_compatible(self):
+#         """
+#         Checks if the underlying RS/CA are compatible
+#         with parameters and concentrations
+#         """
+#         if not isinstance(self.rs, ReactionSystemWithConcentrationsParam):
+#             return False
+#         if not isinstance(self.ca, ContextAutomatonWithConcentrations):
+#             return False
+#         return True
+#
+#     def is_with_concentrations(self):
+#         """
+#         Checks if the underlying RS and CA provide
+#         concentration levels
+#         """
+#         if not isinstance(self.rs, ReactionSystemWithConcentrations):
+#             return False
+#         if not isinstance(self.ca, ContextAutomatonWithConcentrations):
+#             return False
+#         return True
+#
+#     def sanity_check(self):
+#         pass
+#
+#     def get_ordinary_reaction_system_with_automaton(self):
+#
+#         if not self.is_with_concentrations():
+#             raise RuntimeError("Not RS/CA with concentrations")
+#
+#         ors = self.rs.get_reaction_system()
+#         oca = self.ca.get_automaton_with_flat_contexts(ors)
+#
+#         return ReactionSystemWithAutomaton(ors, oca)
+#
+
+# class ReactionSystemWithConcentrationWithAutomaton(ReactionSystemWithAutomaton):
+#
+#     def __init__(self, reaction_system, context_automaton):
+#         self.rs = reaction_system
+#         self.ca = context_automaton
+#
+#     def show(self, soft=False):
+#         self.rs.show(soft)
+#         self.ca.show()

rs_examples.py

@@ -187,3 +187,4 @@ def heat_shock_response(print_system=True,verify_rsc=True):
 def state_translate_rsc2rs(p):
     return [e[0] + "#" + str(e[1]) for e in p]
         
+# EOF

rs_testing.py

@@ -2,16 +2,336 @@ from rs import *
 from smt import *
 import rs_examples 
 from logics import *
+from rsltl_shortcuts import *
+
 import sys
 import resource
 
-def run_tests():
+def run_tests(cmd_args):
     
     # test_extended_automaton()
     # process()
     # heat_shock_response()
     # scalable_chain(print_system=True)
-    example44()
+    # example44()
+    # example44_param()
+    # heat_shock_response_param(cmd_args)
+    # simple_param(cmd_args)
+
+    gene_expression(cmd_args)
+
+def gene_expression(cmd_args):
+    """
+    Simple gene expression example
+    """
+    r = ReactionSystemWithConcentrationsParam()
+    
+    r.add_bg_set_entity(("x", 1))
+    r.add_bg_set_entity(("xp", 1))
+    r.add_bg_set_entity(("X", 1))
+
+    r.add_bg_set_entity(("y", 1))
+    r.add_bg_set_entity(("yp", 1))
+    r.add_bg_set_entity(("Y", 1))
+    
+    r.add_bg_set_entity(("z", 1))
+    r.add_bg_set_entity(("zp", 1))
+    r.add_bg_set_entity(("Z", 1))
+    
+    r.add_bg_set_entity(("h", 1))
+    r.add_bg_set_entity(("Q", 1))
+    r.add_bg_set_entity(("U", 1))
+
+    all_entities = set(r.background_set)
+    
+    r.add_reaction([("x",1)],[("h",1)],[("x",1)])
+
+    # r.add_reaction([("x",1)],[("h",1)],[("xp",1)])
+    param_p1 = r.get_param("P1")
+    r.add_reaction(param_p1,[("h",1)],[("xp",1)])
+    
+    # r.add_reaction([("x",1),("xp",1)],[("h",1)],[("X",1)])
+    param_p2_r = r.get_param("P2_r")
+    param_p2_i = r.get_param("P2_i")
+    param_p2_p = r.get_param("P2_p")
+    r.add_reaction(param_p2_r,param_p2_i,param_p2_p)
+
+    r.add_reaction([("y",1)],[("h",1)],[("y",1)])
+    r.add_reaction([("y",1)],[("Q",1)],[("yp",1)])
+    r.add_reaction([("y",1),("yp",1)],[("h",1)],[("Y",1)])
+    r.add_reaction([("z",1)],[("h",1)],[("z",1)])
+    r.add_reaction([("z",1)], [("X",1)], [("zp",1)])
+    r.add_reaction([("z",1),("zp",1)],[("h",1)],[("Z",1)])
+    r.add_reaction([("U",1),("X",1)],[("h",1)],[("Q",1)])
+
+    c = ContextAutomatonWithConcentrations(r)
+    c.add_init_state("0")
+    c.add_state("1")
+    
+    # the experiments starts with adding x and y:
+    c.add_transition("0", [("x", 1),("y", 1)], "1")
+    
+    # for all the remaining steps we have empty context sequences
+    c.add_transition("1", [], "1")
+
+    rc = ReactionSystemWithAutomaton(r, c)
+    rc.show()
+    
+    f_x1 = ltl_G(True, ltl_Implies(
+        exact_state(["x"], all_entities),
+        ltl_X(bag_And(bag_Not("Y"), bag_Not("Z"), bag_Not("h")), exact_state(["x", "xp"], all_entities))))
+    f_x2 = ltl_G(True, ltl_Implies(
+        exact_state(["x", "xp"], all_entities),
+        ltl_X(bag_Not("h"), exact_state("X", all_entities))))
+    
+    reach_xp = ltl_F(True, "xp")
+    reach_X = ltl_F(True, "X")
+        
+    smt_rsc = SmtCheckerRSCParam(rc, optimise=cmd_args.optimise)
+    
+    smt_rsc.check_rsltl(formulae_list=[f_x1, f_x2, reach_xp, reach_X])
+
+    # smt_rsc.check_rsltl(formula=f_x1, print_witness=True)
+    
+    
+def trivial_param():
+
+    r = ReactionSystemWithConcentrationsParam()
+    r.add_bg_set_entity(("x", 3))
+    r.add_bg_set_entity(("c", 3))
+    r.add_bg_set_entity(("final", 1))
+    
+    param_p1 = r.get_param("P1")
+
+    r.add_reaction(param_p1, [("c", 1)], [("final", 1)])
+    # r.add_reaction([("x", 1)], [("c", 1)], [("final", 1)])
+    
+    c = ContextAutomatonWithConcentrations(r)
+    c.add_init_state("0")
+    c.add_state("1")
+    c.add_transition("0", [("x", 1)], "1")
+    c.add_transition("1", [], "1")
+
+    rc = ReactionSystemWithAutomaton(r, c)
+    rc.show()
+    smt_rsc = SmtCheckerRSCParam(rc)
+
+    f1 = ltl_F(True, bag_entity("final") >= 1)
+
+    # f1 = Formula_rsLTL.f_F(
+    #     BagDescription.f_TRUE(),
+    #     BagDescription.f_entity("final") >= 1)
+    
+    #
+    # WARNING: depth limit is set
+    #
+    smt_rsc.check_rsltl(formula=f1, max_level=10, print_witness=True)
+
+
+def simple_param(cmd_args):
+
+    r = ReactionSystemWithConcentrationsParam()
+    r.add_bg_set_entity(("x", 3))
+    r.add_bg_set_entity(("y", 3))
+    r.add_bg_set_entity(("c", 3))
+    r.add_bg_set_entity(("z", 3))
+    r.add_bg_set_entity(("final", 1))
+    
+    param_p1 = r.get_param("P1")
+
+    r.add_reaction([("x", 1)], [("c", 1)], [("y", 2)])
+    # r.add_reaction([("y", 1)], [("c", 1)], [("z", 1)])
+    r.add_reaction([("y", 1)], [("c", 1)], r.get_param("P2"))
+    r.add_reaction(param_p1, [("x", 1)], [("final", 1)])
+    # r.add_reaction([("z", 1)], [("x", 1)], [("final", 1)])
+
+    c = ContextAutomatonWithConcentrations(r)
+    c.add_init_state("0")
+    c.add_state("1")
+    c.add_transition("0", [("x", 1)], "1")
+    c.add_transition("1", [], "1")
+
+    rc = ReactionSystemWithAutomaton(r, c)
+    rc.show()
+    
+    smt_rsc = SmtCheckerRSCParam(rc, optimise=cmd_args.optimise)
+
+    f1 = Formula_rsLTL.f_F(
+        BagDescription.f_TRUE(),
+        BagDescription.f_entity("final") >= 1)
+    
+    #
+    # WARNING: depth limit is set
+    #
+    smt_rsc.check_rsltl(formula=f1, max_level=10, print_witness=True, cont_if_sat=False)
+
+
+def heat_shock_response_param(cmd_args, print_system=True):
+
+    stress_temp = 42
+    max_temp = 50
+    
+    r = ReactionSystemWithConcentrationsParam()
+    r.add_bg_set_entity(("hsp",1))
+    r.add_bg_set_entity(("hsf",1))
+    r.add_bg_set_entity(("hsf2",1))
+    r.add_bg_set_entity(("hsf3",1))
+    r.add_bg_set_entity(("hse",1))
+    r.add_bg_set_entity(("mfp",1))
+    r.add_bg_set_entity(("prot",1))
+    r.add_bg_set_entity(("hsf3:hse",1))
+    r.add_bg_set_entity(("hsp:mfp",1))
+    r.add_bg_set_entity(("hsp:hsf",1))
+    r.add_bg_set_entity(("stress",1))
+    r.add_bg_set_entity(("no_stress",1))
+    
+    param_p1 = r.get_param("P1")
+    
+    r.add_reaction([("hsf",1)],                         [("hsp",1)],        [("hsf3",1)])
+    r.add_reaction([("hsf",1),("hsp",1),("mfp",1)],     [],                 [("hsf3",1)])
+    r.add_reaction([("hsf3",1)],                        [("hse",1),("hsp",1)],[("hsf",1)])
+    r.add_reaction([("hsf3",1),("hsp",1),("mfp",1)],    [("hse",1)],        [("hsf",1)])
+    r.add_reaction([("hsf3",1),("hse",1)],              [("hsp",1)],        [("hsf3:hse",1)])
+    r.add_reaction([("hsp",1),("hsf3",1),("mfp",1),("hse",1)],[],           [("hsf3:hse",1)])
+    r.add_reaction([("hse",1)],                         [("hsf3",1)],       [("hse",1)])
+    r.add_reaction([("hsp",1),("hsf3",1),("hse",1)],    [("mfp",1)],        [("hse",1)])
+    r.add_reaction([("hsf3:hse",1)],                    [("hsp",1)],        [("hsp",1),("hsf3:hse",1)])
+    r.add_reaction([("hsp",1),("mfp",1),("hsf3:hse",1)],[],                 [("hsp",1),("hsf3:hse",1)])
+    r.add_reaction([("hsf",1),("hsp",1)],               [("mfp",1)],        [("hsp:hsf",1)])
+    r.add_reaction([("hsp:hsf",1),("stress",1)],[],                 [("hsf",1),("hsp",1)])
+    r.add_reaction([("hsp:hsf",1)],                     [("stress",1)],[("hsp:hsf",1)])
+    r.add_reaction([("hsp",1),("hsf3:hse",1)],          [("mfp",1)],        [("hse",1),("hsp:hsf",1)])
+    r.add_reaction([("stress",1),("prot",1)],   [],                     [("mfp",1),("prot",1)])
+    r.add_reaction([("prot",1)],                        [("stress",1)], [("prot",1)])
+    r.add_reaction([("hsp",1),("mfp",1)],               [],                 [("hsp:mfp",1)])
+    r.add_reaction([("mfp",1)],                         [("hsp",1)],        [("mfp",1)])
+    r.add_reaction([("hsp:mfp",1)],                     [],                 [("hsp",1),("prot",1)])
+
+    # r.add_reaction_inc("temp", "heat", [("temp",1)],[("temp",max_temp)])
+    # r.add_reaction_dec("temp", "cool", [("temp",1)],[])
+
+    # r.add_permanency("temp",[("heat",1),("cool",1)])
+
+    c = ContextAutomatonWithConcentrations(r)
+    c.add_init_state("0")
+    c.add_state("1")
+    c.add_transition("0", [("hsf",1),("prot",1),("hse",1),("no_stress",1)], "1")
+    c.add_transition("1", [("stress",1)], "1")
+    c.add_transition("1", [("no_stress",1)], "1")
+    c.add_transition("1", [], "1")
+
+    rc = ReactionSystemWithAutomaton(r,c)
+    
+    if print_system:
+        rc.show()
+    
+    # prop_req = [("hsp:hsf",1),("hse",1),("prot",1)]
+    # prop_block = [("temp",stress_temp)]
+    # prop_req   = [ ("mfp",1) ]
+    # prop_block = [ ]
+    # prop       = (prop_req,prop_block)
+    # rs_prop    = (state_translate_rsc2rs(prop_req),state_translate_rsc2rs(prop_block))
+    #
+    
+    # f_reach_mfp = Formula_rsLTL.f_F(BagDescription.f_TRUE(), (BagDescription.f_entity("mfp") > 0) )
+    f_reach_mfp = ltl_F(True, bag_entity("mfp") > 0)
+
+    f_reach_hspmfp = Formula_rsLTL.f_F(BagDescription.f_TRUE(), (BagDescription.f_entity("hsp:mfp") > 0) )
+    
+    smt_rsc = SmtCheckerRSCParam(rc, optimise=cmd_args.optimise)
+    
+    smt_rsc.check_rsltl(formulae_list=[f_reach_mfp, f_reach_hspmfp])
+    
+    
+    # (1) if we keep increasing the temperature the protein will eventually misfold
+    # f_mfp_when_heating = Formula_rsLTL.f_X(BagDescription.f_TRUE(),
+    #         Formula_rsLTL.f_F(BagDescription.f_entity("heat") > 0, (BagDescription.f_entity("mfp") > 0) )
+    #     )
+    # smt_rsc.check_rsltl(formula=f_mfp_when_heating)
+    
+    # (2) when heating, we finally exceed the stress_temp
+    # f_2 = Formula_rsLTL.f_X(BagDescription.f_TRUE(),
+    #         Formula_rsLTL.f_F(BagDescription.f_entity("heat") > 0, (BagDescription.f_entity("temp") > stress_temp))
+    #     )
+    # smt_rsc.check_rsltl(formula=f_2)
+    
+    # smt_rsc.check_reachability(prop,max_level=40)
+
+def example44_param():
+
+    r = ReactionSystemWithConcentrationsParam()
+    r.add_bg_set_entity(("x", 2))
+    r.add_bg_set_entity(("y", 4))
+    r.add_bg_set_entity(("h", 2))
+    r.add_bg_set_entity(("m", 1))
+
+    r.add_reaction([("y", 1), ("x", 1)], [("y", 2), ("h", 1)], [("y", 2)])
+    r.add_reaction([("y", 2), ("x", 2)], [("y", 3), ("h", 1)], [("y", 3)])
+    r.add_reaction([("y", 3), ("h", 1)], [("y", 4), ("h", 2)], [("y", 4)])
+    r.add_reaction([("y", 4), ("h", 1)], [("h", 2)], [("y", 3)])
+    r.add_reaction([("y", 4), ("x", 2)], [("h", 1)], [("y", 2)])
+    r.add_reaction([("m", 1)], [("y", 3)], [("m", 1)])
+
+    c = ContextAutomatonWithConcentrations(r)
+    c.add_init_state("0")
+    c.add_state("1")
+    c.add_transition("0", [("y", 1), ("m", 1), ("x", 1)], "1")
+    c.add_transition("1", [("x", 1)], "1")
+    c.add_transition("1", [("x", 1), ("h", 1)], "1")
+    c.add_transition("1", [("x", 2)], "1")
+    c.add_transition("1", [("x", 2), ("h", 1)], "1")
+    c.add_transition("1", [("h", 1)], "1")
+
+    rc = ReactionSystemWithAutomaton(r, c)
+    rc.show()
+    smt_rsc = SmtCheckerRSCParam(rc)
+
+    # Universal property which seems to be true: (holds also existentially)
+    f1 = Formula_rsLTL.f_G(BagDescription.f_entity("x") > 0,
+                           Formula_rsLTL.f_Implies(
+        (BagDescription.f_entity('y') == 2),
+        Formula_rsLTL.f_X(
+            (BagDescription.f_entity("x") > 1),
+            BagDescription.f_entity("y") >= 3
+        )
+    )
+    )
+
+    # lets see if we can find a counterexample to this property:
+    neg_f1 = Formula_rsLTL.f_F(BagDescription.f_entity("x") > 0,
+                               Formula_rsLTL.f_And(
+        (BagDescription.f_entity('y') == 2),
+        Formula_rsLTL.f_X(
+            (BagDescription.f_entity("x") > 1),
+            BagDescription.f_entity("y") < 3
+        )
+    )
+    )
+
+    # we fix the property f1
+    # this one holds:
+    f2 = Formula_rsLTL.f_G(
+        BagDescription.f_entity("x") > 0, Formula_rsLTL.f_Implies(
+            (BagDescription.f_entity('y') == 2),
+            Formula_rsLTL.f_X(
+                (BagDescription.f_entity("x") > 1) &
+                (BagDescription.f_entity("h") < 1),
+                BagDescription.f_entity("y") >= 3)))
+
+    # neg_f1 = Formula_rsLTL.f_X(BagDescription.f_TRUE(), Formula_rsLTL.f_F(BagDescription.f_entity("x") > 0,
+    #     Formula_rsLTL.f_And(
+    #         (BagDescription.f_entity('y') > 0),
+    #         Formula_rsLTL.f_X(
+    #             BagDescription.f_entity("x") > 0,
+    #             Formula_rsLTL.f_X(
+    #                 BagDescription.f_entity("x") > 1,
+    #                 BagDescription.f_entity("y") < 3
+    #             )
+    #         )
+    #     )
+    # ))
+    smt_rsc.check_rsltl(formula=neg_f1)
+
 
 def example44():
 
@@ -38,7 +358,6 @@ def example44():
     c.add_transition("1", [("x", 2), ("h", 1)], "1")
     c.add_transition("1", [("h", 1)], "1")
 
-
     rc = ReactionSystemWithAutomaton(r, c)
     rc.show()
     smt_rsc = SmtCheckerRSC(rc)
@@ -67,15 +386,13 @@ def example44():
 
     # we fix the property f1
     # this one holds:
-    f2 = Formula_rsLTL.f_G(BagDescription.f_entity("x") > 0, 
-        Formula_rsLTL.f_Implies(
+    f2 = Formula_rsLTL.f_G(
+        BagDescription.f_entity("x") > 0, Formula_rsLTL.f_Implies(
             (BagDescription.f_entity('y') == 2),
             Formula_rsLTL.f_X(
-                (BagDescription.f_entity("x") > 1) & (BagDescription.f_entity("h") < 1), 
-                BagDescription.f_entity("y") >= 3
-            )
-        )
-    )
+                (BagDescription.f_entity("x") > 1) &
+                (BagDescription.f_entity("h") < 1),
+                BagDescription.f_entity("y") >= 3)))
 
     # neg_f1 = Formula_rsLTL.f_X(BagDescription.f_TRUE(), Formula_rsLTL.f_F(BagDescription.f_entity("x") > 0,
     #     Formula_rsLTL.f_And(
@@ -91,6 +408,7 @@ def example44():
     # ))
     smt_rsc.check_rsltl(formula=neg_f1)
 
+    
 def heat_shock_response(print_system=True):
 
     stress_temp = 42

rsltl_shortcuts.py

@@ -0,0 +1,75 @@
+from logics import *
+
+##### SHORTCUTS
+
+def bag_True():
+    return BagDescription.f_TRUE()
+
+def bag_entity(name):
+    return BagDescription.f_entity(name)
+
+def get_bag_if_str(arg):
+    if isinstance(arg, str):
+        return bag_entity(arg) > 0
+    else:
+        return arg
+
+def bag_Not(a0):
+    a0 = get_bag_if_str(a0)
+    return BagDescription.f_Not(a0)
+
+def bag_And(*args):
+    assert len(args) > 1
+    last = get_bag_if_str(args[0])
+    for arg in args[1:]:
+        last = BagDescription.f_And(
+            last, get_bag_if_str(arg))
+    return last
+    
+def exact_state(contained_entities, all_entities):
+    expr = []
+    for ent_str in all_entities:
+        ent = bag_entity(ent_str)
+        if ent_str in contained_entities:
+            expr.append(ent > 0)
+        else:
+            expr.append(ent == 0)
+
+    if len(expr) > 0:
+        
+        last = expr[0]
+        for e in expr[1:]:
+            last = BagDescription.f_And(last, e)
+        return last
+        
+    else:
+        assert False
+
+def ltl_F(ctx_arg, a0):
+    a0 = get_bag_if_str(a0)
+    return Formula_rsLTL.f_F(ctx_arg, a0)
+
+def ltl_G(ctx_arg, a0):
+    a0 = get_bag_if_str(a0)
+    return Formula_rsLTL.f_G(ctx_arg, a0)
+
+def ltl_X(ctx_arg, a0):
+    a0 = get_bag_if_str(a0)
+    return Formula_rsLTL.f_X(ctx_arg, a0)
+
+def ltl_And(*args):
+    assert len(args) > 1
+    last = get_bag_if_str(args[0])
+    for arg in args[1:]:
+        last = Formula_rsLTL.f_And(last, get_bag_if_str(arg))
+    return last
+
+def ltl_Not(a0):
+    a0 = get_bag_if_str(a0)
+    return Formula_rsLTL.f_Not(a0)
+    
+def ltl_Implies(a0, a1):
+    a0 = get_bag_if_str(a0)
+    a1 = get_bag_if_str(a1)
+    return Formula_rsLTL.f_Implies(a0, a1)
+

rssmt.py

@@ -6,6 +6,8 @@
 
 """
 
+import argparse
+
 from rs import *
 from smt import *
 import sys
@@ -14,10 +16,11 @@ import rs_testing
 
 from colour import *
 
-import resource
-
 profiling = False
 
+if profiling:
+    import resource
+
 ##################################################################
 
 version = "2017/04/08/00"
@@ -25,11 +28,12 @@ rsmc_banner = """
 Reaction Systems SMT-Based Model Checking
 
 Version:  """ + version + """
-Author:   Artur Męski <meski@ipipan.waw.pl> / <artur.meski@ncl.ac.uk>
+Author:   Artur Meski <meski@ipipan.waw.pl> / <artur.meski@ncl.ac.uk>
 """
 
 ##################################################################
 
+
 def print_banner():
     print()
     for line in rsmc_banner.split("\n"):
@@ -38,11 +42,21 @@ def print_banner():
 
 ##################################################################
 
+
 def main():
     """Main function"""
 
+    parser = argparse.ArgumentParser()
+
+    parser.add_argument("-v", "--verbose", 
+        help="turn verbosity on", action="store_true")
+    parser.add_argument("-o", "--optimise", 
+        help="minimise the parametric computation result", action="store_true")
+
+    args = parser.parse_args()
+
     print_banner()
-    rs_testing.run_tests()
+    rs_testing.run_tests(args)
 
 ##################################################################
 

smt/__init__.py

@@ -1,5 +1,6 @@
 #from smt.smt_checker import SmtChecker
 from smt.smt_checker_rs import SmtCheckerRS
 from smt.smt_checker_rsc import SmtCheckerRSC
+from smt.smt_checker_rsc_param import SmtCheckerRSCParam
 from smt.smt_checker_rs_na import SmtCheckerRSNA
 from smt.smt_checker_distrib_rs import SmtCheckerDistribRS

smt/smt_checker_rsc_param.py

@@ -0,0 +1,961 @@
+"""
+SMT-based Model Checking Module for RS with Concentrations and Context Automaton
+"""
+
+from z3 import *
+from time import time
+from sys import stdout
+from itertools import chain
+import resource
+from colour import *
+
+from logics import rsLTL_Encoder
+
+from rs.reaction_system_with_concentrations_param import ParameterObj, is_param
+
+# def simplify(x):
+#     return x
+
+def z3_max(a, b):
+    return If(a > b, a, b)
+
+class SmtCheckerRSCParam(object):
+
+    def __init__(self, rsca, optimise=False):
+    
+        rsca.sanity_check()
+
+        if not rsca.is_concentr_and_param_compatible():
+            raise RuntimeError("RS and CA with concentrations (and parameters) expected")
+        
+        self.rs = rsca.rs
+        self.ca = rsca.ca
+        
+        self.optimise = optimise
+
+        self.initialise()
+    
+    def initialise(self):
+        """Initialises all the variables used by the checker"""
+
+        ### "Currently" used variables (loaded from self.path_v...)
+        self.v = None
+        self.v_ctx = None
+        self.ca_state = None
+        
+        # intermediate products:
+        self.v_improd = None
+        self.v_improd_for_entities = None
+        
+        ### Per-path variables
+        self.path_v = dict()
+        self.path_v_ctx = dict()
+        self.path_ca_state = dict()
+        
+        # intermediate products:
+        self.path_v_improd = dict()
+        self.path_v_improd_for_entities = dict()
+        
+        # parameters:
+        self.v_param = dict()
+
+        self.next_level_to_encode = 0
+
+        # this is probably not needed anymore:
+        # self.producible_entities = self.rs.get_producible_entities()
+        # self.improducible_entities = set(self.rs.get_state_ids(self.rs.background_set)) - self.producible_entities
+        
+        # WARNING: improd vs. improducible
+        # there is some confusion related to the variables naming:
+        # improd - intermediate products
+        # improducible - entities that are never produces (there is no reaction that produces that entity)
+        
+        # TODO: number of loops == number of paths
+        self.loop_position = None
+        self.path_loop_position = dict()
+        
+        if self.optimise:
+            self.solver = Optimize()
+        else:
+            self.solver = Solver() #For("QF_FD")
+
+        self.verification_time = None
+                
+        self.prepare_param_variables()
+
+    def reset(self):
+        """Reinitialises the state of the checker"""
+        
+        self.initialise()
+        
+    def prepare_all_variables(self, num_of_paths):
+
+        for path_idx in range(num_of_paths):
+            self.prepare_all_path_variables(path_idx)
+        self.next_level_to_encode += 1
+        
+    def prepare_all_path_variables(self, path_idx):
+        """Prepares the variables for a given path index"""
+
+        print_info("Preparing variables for path={:d} (level={:d})".format(path_idx, self.next_level_to_encode))
+        
+        self.prepare_state_variables(path_idx)
+        self.prepare_context_variables(path_idx)
+        self.prepare_intermediate_product_variables(path_idx)
+        self.prepare_loop_position_variables(path_idx)
+            
+    def prepare_loop_position_variables(self, path_idx):
+        """Prepares the variables for loop positions"""
+        
+        self.path_loop_position[path_idx] = Int("p{:d}_loop_pos".format(path_idx))
+        
+    def prepare_context_variables(self, path_idx):
+        """Prepares all the context variables"""
+
+        level = self.next_level_to_encode
+
+        self.path_v_ctx.setdefault(path_idx, [])
+        assert len(self.path_v_ctx[path_idx]) == level
+        
+        variables = []
+        for entity in self.rs.background_set:
+            new_var = Int("p{:d}C{:d}_{:s}".format(path_idx, level, entity))
+            variables.append(new_var)
+
+        self.path_v_ctx[path_idx].append(variables)
+
+    def prepare_state_variables(self, path_idx):
+        """Prepares all the state variables"""
+
+        level = self.next_level_to_encode
+        
+        # RS vars
+        self.path_v.setdefault(path_idx, [])
+        assert len(self.path_v[path_idx]) == level
+
+        variables = []
+        for entity in self.rs.background_set:
+            new_var = Int("p{:d}L{:d}_{:s}".format(path_idx, level, entity))
+            variables.append(new_var)
+        self.path_v[path_idx].append(variables)
+        
+        # Context automaton states:
+        self.path_ca_state.setdefault(path_idx, [])
+        assert len(self.path_ca_state[path_idx]) == level
+        
+        ca_state_var = Int("p{:d}CA{:d}_state".format(path_idx, level))
+        self.path_ca_state[path_idx].append(ca_state_var)
+
+    def prepare_intermediate_product_variables(self, path_idx):
+        """
+        Prepares the intermediate product variables
+        carrying the individual concentration levels produced
+        by the reactions.
+
+        These variables are used later on to encode the final
+        concentration levels for all the entities
+        """
+
+        level = self.next_level_to_encode
+
+        self.path_v_improd.setdefault(path_idx, [])
+        self.path_v_improd_for_entities.setdefault(path_idx, [])
+        
+        assert len(self.path_v_improd[path_idx]) == level
+        assert len(self.path_v_improd_for_entities[path_idx]) == level
+
+        if level < 1:
+            #
+            # If we are at level==0, we add a dummy "level"
+            # to match the indices of of the successors
+            # which are always at level+1.
+            #
+            self.path_v_improd[path_idx].append(None)
+            self.path_v_improd_for_entities[path_idx].append(None)
+        else:
+
+            reactions_dict = dict()
+            number_of_reactions = len(self.rs.reactions)
+
+            all_entities_dict = dict()
+
+            for reaction in self.rs.reactions:
+                *_, products = reaction
+
+                reaction_id = self.rs.reactions.index(reaction)
+
+                entities_dict = dict()
+                
+                if is_param(products):
+
+                    for entity in self.rs.set_of_bgset_ids:
+                        entity_name = self.rs.get_entity_name(entity)
+                        new_var = Int("p{:d}L{:d}_ImProd_r{:d}_{:s}".format(path_idx, level, reaction_id, entity_name))
+                        entities_dict[entity] = new_var
+
+                        all_entities_dict.setdefault(entity, [])
+                        all_entities_dict[entity].append(new_var)
+                                            
+                else:
+                
+                    for entity, conc in products:
+                        entity_name = self.rs.get_entity_name(entity)
+                        new_var = Int("p{:d}L{:d}_ImProd_r{:d}_{:s}".format(path_idx, level, reaction_id, entity_name))
+                        entities_dict[entity] = new_var
+
+                        all_entities_dict.setdefault(entity, [])
+                        all_entities_dict[entity].append(new_var)
+
+                reactions_dict[reaction_id] = entities_dict
+            
+            self.path_v_improd[path_idx].append(reactions_dict)
+            self.path_v_improd_for_entities[path_idx].append(all_entities_dict)
+        
+    def prepare_param_variables(self):
+        """
+        Prepares variables for parameters
+        
+        A parameter (it's valuation) is a subset of the background set,
+        therefore we need separate variables for each element of the
+        background set.
+        """
+        
+        for param_name in self.rs.parameters.keys():
+
+            # we start collecting bg-related vars for the given param
+            vars_for_param = []
+
+            for entity in self.rs.set_of_bgset_ids:
+                new_var = Int("Pm{:s}_{:s}".format(param_name, self.rs.get_entity_name(entity)))
+                vars_for_param.append(new_var)
+
+            self.v_param[param_name] = vars_for_param
+            
+    def load_varset_for_path(self, path_idx):
+        """
+        Loads the the variables for the path with path_idx
+        """
+        self.v = self.path_v[path_idx]
+        self.v_ctx = self.path_v_ctx[path_idx]
+        self.ca_state = self.path_ca_state[path_idx]        
+        self.v_improd = self.path_v_improd[path_idx]
+        self.v_improd_for_entities = self.path_v_improd_for_entities[path_idx]
+        self.loop_position = self.path_loop_position[path_idx]
+
+    def enc_param_concentration_levels_assertion(self):
+        """
+        Assertions for the parameter variables
+        """
+        
+        if len(self.v_param) == 0:
+            return True
+        
+        enc_param_gz = True
+        enc_non_empty = True
+        
+        for param_vars in self.v_param.values():
+            
+            enc_param_at_least_one = False
+            for pvar in param_vars:
+            
+                # TODO: fixed upper limit: 100 (have a per-param setting for that)
+                enc_param_gz = simplify(And(enc_param_gz, pvar >= 0, pvar < 100))
+                enc_param_at_least_one = simplify(Or(enc_param_at_least_one, pvar > 0))
+            
+            enc_non_empty = simplify(And(enc_non_empty, enc_param_at_least_one))
+            
+        return simplify(And(enc_param_gz, enc_non_empty))
+    
+    def assert_param_optimisation(self):
+
+        for param_vars in self.v_param.values():
+            for pvar in param_vars:
+                self.solver.add_soft(pvar < 1)
+                                        
+    def enc_concentration_levels_assertion(self, level, path_idx):
+        """
+        Encodes assertions that (some) variables need to be >=0
+        
+        We do not need to actually control all the variables,
+        only those that can possibly go below 0.
+        """
+
+        print_info("Concentration level assertions for path={:d} (level={:d})".format(path_idx, level))
+        
+        enc_gz = True
+
+        for e_i in self.rs.set_of_bgset_ids:
+            var = self.path_v[path_idx][level][e_i]
+            var_ctx = self.path_v_ctx[path_idx][level][e_i]
+            e_max = self.rs.get_max_concentration_level(e_i)
+            enc_gz = simplify(And(enc_gz, var >= 0, var_ctx >= 0, var <= e_max, var_ctx <= e_max))
+            
+            vars_per_reaction = self.path_v_improd_for_entities[path_idx][level + 1]
+            if e_i in vars_per_reaction:
+                for var_improd in vars_per_reaction[e_i]:
+                    enc_gz = simplify(And(enc_gz, var_improd >= 0, var_improd <= e_max))
+        
+        return enc_gz
+ 
+    def enc_init_state(self, level, path_idx):
+        """Encodes the initial state at the given level"""
+
+        rs_init_state_enc = True
+
+        for v in self.path_v[path_idx][level]:
+            # the initial concentration levels are zeroed
+            rs_init_state_enc = simplify(And(rs_init_state_enc, v == 0))
+
+        ca_init_state_enc = self.path_ca_state[path_idx][level] == self.ca.get_init_state_id()
+
+        init_state_enc = simplify(And(rs_init_state_enc, ca_init_state_enc))
+
+        return init_state_enc
+
+    def enc_transition_relation(self, level, path_idx):
+        return simplify(
+            And(self.enc_rs_trans(level, path_idx),
+                self.enc_automaton_trans(level, path_idx)))
+
+    def enc_single_reaction(self, level, path_idx, reaction):
+        """
+        Encodes a single reaction
+
+        For encoding the products we use intermediate variables:
+
+            * each reaction has its own product variables,
+
+            * those are meant to be used to compute the MAX concentration
+
+        """
+
+        reactants, inhibitors, products = reaction
+
+        # we need reaction_id to find the intermediate product variable
+        reaction_id = self.rs.reactions.index(reaction)
+
+        # ** REACTANTS *******************************************
+        enc_reactants = True
+        if is_param(reactants):
+            param_name = reactants.name
+            for entity in self.rs.set_of_bgset_ids:
+                enc_reactants = And(enc_reactants, Or(
+                    self.v_param[param_name][entity] == 0,
+                    self.path_v[path_idx][level][entity] >= self.v_param[param_name][entity], 
+                    self.path_v_ctx[path_idx][level][entity] >= self.v_param[param_name][entity]))
+        else:
+            for entity, conc in reactants:
+                enc_reactants = And(enc_reactants, Or(
+                    self.path_v[path_idx][level][entity] >= conc, 
+                    self.path_v_ctx[path_idx][level][entity] >= conc))
+
+        # ** INHIBITORS ******************************************
+        enc_inhibitors = True
+        if is_param(inhibitors):
+            param_name = inhibitors.name
+            for entity in self.rs.set_of_bgset_ids:
+                enc_inhibitors = And(enc_inhibitors, Or(
+                    self.v_param[param_name][entity] == 0,
+                    And(
+                        self.path_v[path_idx][level][entity] < self.v_param[param_name][entity], 
+                        self.path_v_ctx[path_idx][level][entity] < self.v_param[param_name][entity])))
+        else:
+            for entity, conc in inhibitors:
+                enc_inhibitors = And(enc_inhibitors, And(
+                    self.path_v[path_idx][level][entity] < conc, 
+                    self.path_v_ctx[path_idx][level][entity] < conc))
+
+        # ** PRODUCTS *******************************************
+        enc_products = True
+        if is_param(products):
+            param_name = products.name
+            for entity in self.rs.set_of_bgset_ids:
+                enc_products = simplify(And(enc_products,
+                    self.path_v_improd[path_idx][level + 1][reaction_id][entity] == self.v_param[param_name][entity]))
+        else:
+            for entity, conc in products:
+                enc_products = simplify(And(enc_products, 
+                    self.path_v_improd[path_idx][level + 1][reaction_id][entity] == conc))
+
+        # Nothing is produced (when the reaction is disabled)
+        enc_no_prod = True
+        if is_param(products):
+            for entity in self.rs.set_of_bgset_ids:
+                enc_no_prod = And(enc_no_prod,
+                    self.path_v_improd[path_idx][level + 1][reaction_id][entity] == 0)
+        else:
+            for entity, _ in products:
+                enc_no_prod = simplify(And(enc_no_prod, 
+                    self.path_v_improd[path_idx][level + 1][reaction_id][entity] == 0))
+        
+        #
+        # (R and I) iff P
+        #
+        enc_enabled = And(enc_reactants, enc_inhibitors) == enc_products
+        
+        #
+        # ~(R and I) iff P_zero
+        #
+        enc_not_enabled = Not(And(enc_reactants, enc_inhibitors)) == enc_no_prod
+        
+        enc_reaction = And(enc_enabled, enc_not_enabled)
+
+        return enc_reaction
+
+    def enc_general_reaction_enabledness(self, level, path_idx):
+        """
+        General enabledness condition for reactions
+        
+        The necessary condition for a reaction to be enabled is
+        that the state is not empty, i.e., at least one entity
+        is present in the current state.
+        
+        This condition must be used when there are parametric reactions
+        because parameters could have all the entities set to zero and
+        that immediately allows for all the conditions on the reactants
+        to be fulfilled: (entity <= param) -> (0 <= 0)
+        """
+        
+        enc_cond = False
+        for entity in self.rs.set_of_bgset_ids:
+            enc_cond = simplify(Or(enc_cond, self.path_v[path_idx][level][entity] > 0, self.path_v_ctx[path_idx][level][entity] > 0))
+
+        return enc_cond
+
+    def enc_rs_trans(self, level, path_idx):
+        """Encodes the transition relation"""
+
+        #
+        # IMPORTANT NOTE
+        #
+        # We need to make sure we do something about the UNUSED ENTITIES
+        # that is, those that are never produced.
+        #
+        # They should have concentration levels set to 0.
+        #
+        # That needs to happen automatically (in the MAX encoding) -- for the parametric 
+        # case it makes no sense to identify the entities that are never produced, unless 
+        # we have no parameters as products (special case, so that could be an 
+        # optimisation)
+        #
+
+        enc_trans = True
+
+        for reaction in self.rs.reactions:
+            enc_reaction = self.enc_single_reaction(level, path_idx, reaction)
+            enc_trans = simplify(And(enc_trans, enc_reaction))
+
+        # Next we encode the MAX concentration values: 
+        # we collect those from the intermediate product variables
+
+        enc_max_prod = True
+
+        # Save all the intermediate product variables for a given level:
+        #
+        # - Intermediate products of (level+1) correspond to the next level
+        #
+        #   {reactants & inhibitors}[level] 
+        #       => 
+        #   {improd}[level+1] 
+        #       => 
+        #   {products}[level+1]
+        #
+        current_v_improd_for_entities = self.path_v_improd_for_entities[path_idx][level + 1]
+        for entity in self.rs.set_of_bgset_ids:
+            per_reaction_vars = current_v_improd_for_entities.get(entity, [])
+            enc_max_prod = simplify(
+                And(enc_max_prod, self.path_v[path_idx][level + 1][entity] == self.enc_max(per_reaction_vars)))
+        
+        # make sure at least one entity is >0
+        enc_general_cond = self.enc_general_reaction_enabledness(level, path_idx)
+
+        enc_trans_with_max = simplify(And(enc_general_cond, enc_max_prod, enc_trans))
+
+        # print(enc_trans_with_max)
+
+        return enc_trans_with_max
+
+    def enc_max(self, elements):
+
+        enc = None
+        
+        if elements == []:
+            enc = 0
+
+        elif len(elements) == 1:
+            enc = z3_max(0, elements[0])
+
+        elif len(elements) > 1:
+
+            enc = 0
+            for i in range(len(elements) - 1):
+                enc = z3_max(enc, z3_max(elements[i], elements[i + 1]))
+
+        return enc
+    
+    def enc_automaton_trans(self, level, path_idx):
+        """Encodes the transition relation for the context automaton"""
+        
+        enc_trans = False
+    
+        for src,ctx,dst in self.ca.transitions:
+                src_enc = self.path_ca_state[path_idx][level] == src
+                dst_enc = self.path_ca_state[path_idx][level+1] == dst
+                
+                all_ent = set(range(len(self.rs.background_set)))
+                
+                incl_ctx = set([e for e,c in ctx])
+                excl_ctx = all_ent - incl_ctx
+                
+                ctx_enc = True
+                
+                for e,c in ctx:
+                    ctx_enc = simplify(And(ctx_enc, self.path_v_ctx[path_idx][level][e] == c))
+                
+                for e in excl_ctx:
+                    ctx_enc = simplify(And(ctx_enc, self.path_v_ctx[path_idx][level][e] == 0))
+                
+                cur_trans = simplify(And(src_enc, ctx_enc, dst_enc))
+                enc_trans = simplify(Or(enc_trans, cur_trans))
+        
+        return enc_trans
+
+    def enc_exact_state(self, level, state):
+        """Encodes the state at the given level with the exact concentration values"""
+
+        raise RuntimeError("Should not be used with RSC")
+
+    def enc_min_state(self, level, state):
+        """Encodes the state at the given level with the minimal required concentration levels"""
+
+        enc = True
+        for ent,conc in state:
+            e_id = self.rs.get_entity_id(ent)
+            enc = And(enc, self.v[level][e_id] >= conc)
+
+        # state_ids = self.rs.get_state_ids(state)
+        #
+        # for entity in state_ids:
+        #     enc = And(enc, self.v[level][entity])
+
+        return simplify(enc)    
+
+    def enc_state_with_blocking(self, level, prop):
+        """Encodes the state at the given level with blocking certain concentrations"""
+
+        required,blocked = prop
+
+        enc = True
+        for ent,conc in required:
+                e_id = self.rs.get_entity_id(ent)
+                enc = And(enc, self.v[level][e_id] >= conc)
+
+        for ent,conc in blocked:
+                e_id = self.rs.get_entity_id(ent)
+                enc = And(enc, self.v[level][e_id] < conc)
+
+        return simplify(enc)
+
+    def decode_witness(self, max_level, path_idx, print_model=False):
+        """
+        Decodes the witness
+        
+        Also decodes the parameters
+        """
+
+        m = self.solver.model()
+
+        if print_model:
+            print(m)
+
+        for level in range(max_level + 1):
+
+            print("\n{: >70}".format("[ level=" + repr(level) + " ]"))
+
+            print("  State: {", end=""),
+            for var_id in range(len(self.path_v[path_idx][level])):
+                var_rep = repr(m[self.path_v[path_idx][level][var_id]])
+                if not var_rep.isdigit():
+                    raise RuntimeError(
+                        "unexpected: representation is not a positive integer")
+                if int(var_rep) > 0:
+                    print(
+                        " " + self.rs.get_entity_name(var_id) + "=" + var_rep,
+                        end="")
+                # print(" " + repr(m[self.v[level][var_id]]), end="")
+            print(" }")
+
+            if level != max_level:
+                print("  Context set: ", end="")
+                print("{", end="")
+                for var_id in range(len(self.path_v[path_idx][level])):
+                    var_rep = repr(m[self.path_v_ctx[path_idx][level][var_id]])
+                    if not var_rep.isdigit():
+                        raise RuntimeError(
+                            "unexpected: representation is not a positive integer")
+                    if int(var_rep) > 0:
+                        print(
+                            " " + self.rs.get_entity_name(var_id) + "=" + var_rep,
+                            end="")
+                print(" }")
+        
+        # Parameters    
+        
+        print("\n\n  Parameters:\n")
+        for param_name in self.rs.parameters.keys():
+            print("{: >6}: ".format(param_name), end="")
+            print("{", end="")
+            
+            params = self.v_param[param_name]
+            
+            for entity in self.rs.set_of_bgset_ids:
+                var_rep = repr(m[params[entity]])
+                if not var_rep.isdigit():
+                    raise RuntimeError(
+                        "unexpected: representation is not a positive integer")
+                if int(var_rep) > 0:
+                    print(
+                        " " + str(self.rs.get_entity_name(entity)) + "=" + str(var_rep),
+                        end="")
+            print(" }")
+        
+        print("\n")
+            
+    def check_rsltl(
+            self, formulae_list, 
+            print_witness=True, 
+            print_time=True, print_mem=True,
+            max_level=None, cont_if_sat=False):
+        """
+        Bounded Model Checking for rsLTL properties
+            
+            * print_witness   -- prints the decoded witness
+            * print_time      -- prints the time consumed
+            * print_mem       -- prints the memory consumed
+            * max_level       -- if not None, the methods 
+                                 stops at the specified level
+            * cont_if_sat     -- if True, then the method
+                                 continues up until max_level is
+                                 reached (even if sat found)
+        """
+            
+        if not isinstance(formulae_list, (list, tuple)):
+            print_error("Expected a list of formulae")
+
+        self.reset()
+
+        num_of_paths = len(formulae_list)
+        
+        # formula = formulae[0]
+
+        print_info("Running rsLTL bounded model checking")
+        print_info("Tested formulae:")
+        for form in formulae_list:
+            print_info(" "*4 + str(form))
+
+        if print_time:
+            start = resource.getrusage(resource.RUSAGE_SELF).ru_utime
+
+        self.prepare_all_variables(num_of_paths)        
+        
+        self.load_varset_for_path(0)
+        
+        # initial states for all the paths
+        initial_states = []
+        for path_idx in range(num_of_paths):
+            initial_states.append(self.enc_init_state(0, path_idx))
+        self.solver_add(initial_states)
+        
+        self.current_level = 0
+
+        self.prepare_all_variables(num_of_paths)
+
+        # assertions for all the paths and parameters
+        additional_assertions = []
+        for path_idx in range(num_of_paths):
+            additional_assertions.append(self.enc_concentration_levels_assertion(0, path_idx))
+        additional_assertions.append(self.enc_param_concentration_levels_assertion())
+        self.solver_add(additional_assertions)
+
+        if self.optimise:
+            self.assert_param_optimisation()
+
+        encoder = rsLTL_Encoder(self)
+
+        print_info("Iterating...")
+        
+        while True:
+            
+
+            
+            print(
+                "\n{:-^70}".format("[ Working at level=" + str(self.current_level) + " ]"))
+            # stdout.flush()
+
+            # reachability test:
+            self.solver.push()
+
+            enc_form = []
+            for formula in formulae_list:
+                
+                path_idx = formulae_list.index(formula)
+                
+                print_info("Generating the encoding for {:s} ({:d} of {:d})".format(
+                    str(formula), path_idx+1, len(formulae_list)))
+
+                encoder.load_variables(
+                    var_rs=self.path_v[path_idx],
+                    var_ctx=self.path_v_ctx[path_idx],
+                    var_loop_pos=self.path_loop_position[path_idx])
+                    
+                enc_form.append(encoder.get_encoding(formula, self.current_level))
+                ncalls = encoder.get_ncalls()
+
+                print_info("Cache hits: {:d}, encode calls: {:d} (approx: {:d})".format(
+                        encoder.get_cache_hits(), ncalls[0], ncalls[1]))
+
+                encoder.flush_cache()
+                
+            print_info("Adding the formulae to the solver...")
+            
+            # print (enc_form)
+            
+            self.solver_add(enc_form)
+
+            print_info("Adding the loops encoding...")
+            self.solver_add(self.get_loop_encodings())
+
+            print_info("Testing satisfiability...")
+            result = self.solver.check()
+            if result == sat:
+                print_positive(green_str(
+                        "SAT at level={:d}".format(self.current_level)))
+                # print(self.solver.model())
+                if print_witness:
+                    for formula in formulae_list:
+                        path_idx = formulae_list.index(formula)
+                        print("\n{:=^70}".format("[ WITNESS ]"))
+                        print("\n  Witness for: {:s}".format(str(formula)))
+                        self.decode_witness(self.current_level, path_idx)
+                if not cont_if_sat:
+                    break
+            else:
+                print_info("Unsat")
+                self.solver.pop()
+
+            self.prepare_all_variables(num_of_paths)
+            
+            # assertions for all the paths
+            additional_assertions = []
+            for path_idx in range(num_of_paths):
+                additional_assertions.append(
+                    self.enc_concentration_levels_assertion(self.current_level + 1, path_idx))
+            self.solver_add(additional_assertions)
+
+            print_info("Unrolling the transition relation")
+            for path_idx in range(num_of_paths):
+                self.solver_add(self.enc_transition_relation(self.current_level, path_idx))
+
+            print(
+                "{:->70}".format("[ level=" + str(self.current_level) + " done ]"))
+            self.current_level += 1
+
+            if not max_level is None and self.current_level > max_level:
+                print("Stopping at level=" + str(max_level))
+                break
+
+        if print_time:
+            # stop = time()
+            stop = resource.getrusage(resource.RUSAGE_SELF).ru_utime
+            self.verification_time = stop - start
+            print()
+            print(
+                "\n[i] {: >60}".format(
+                    " Time: " + repr(self.verification_time) + " s"))
+
+        if print_mem:
+            print(
+                "[i] {: >60}".format(
+                    " Memory: " +
+                    repr(
+                        resource.getrusage(resource.RUSAGE_SELF).ru_maxrss /
+                        (1024 * 1024)) + " MB"))
+        
+    def dummy_unroll(self, levels):
+        """Unrolls the variables for testing purposes"""
+
+        self.current_level = -1     
+        for i in range(levels+1):
+            self.prepare_all_variables()
+            self.current_level += 1
+
+        print(C_MARK_INFO + " Dummy Unrolling done.")
+        
+    def state_equality(self, level_A, level_B):
+        """Encodes equality of two states at two different levels"""
+        
+        eq_enc = True
+        
+        for e_i in range(len(self.rs.background_set)):
+            e_i_equality = self.v[level_A][e_i] == self.v[level_B][e_i]
+            eq_enc = simplify(And(eq_enc, e_i_equality))
+
+        eq_enc_ctxaut = self.ca_state[level_A] == self.ca_state[level_B]
+        eq_enc = simplify(And(eq_enc, eq_enc_ctxaut))
+
+        return eq_enc   
+    
+    def get_loop_encodings(self):
+        
+        k = self.current_level
+        loop_var = self.loop_position
+        
+        loop_enc = True
+        
+        """
+        (loop_var == i) means that there is a loop taking back to the state (i-1)
+        
+        Therefore, the encoding starts at 1, not at 0.
+        """
+        
+        for i in range(1,k+1):
+            loop_enc = simplify(And(loop_enc, Implies( loop_var == i, self.state_equality(i-1, k) )))
+        
+        return loop_enc
+    
+    def solver_add(self, expression):
+        """
+        This is a solver.add() wrapper
+        """
+    
+        if expression == True:
+            return
+    
+        if expression == False:
+            raise RuntimeError("Trying to assert False.")
+        
+        self.solver.add(expression)
+
+    def check_reachability(self, state, print_witness=True, 
+            print_time=True, print_mem=True, max_level=1000):
+        """Main testing function"""
+
+        self.reset()
+
+        if print_time:
+            # start = time()
+            start = resource.getrusage(resource.RUSAGE_SELF).ru_utime
+
+        self.prepare_all_variables()
+        self.solver_add(self.enc_init_state(0))
+        self.current_level = 0
+
+        self.prepare_all_variables()
+        
+        self.solver_add(self.enc_concentration_levels_assertion(0))
+        
+        while True:
+            self.prepare_all_variables()
+            self.solver_add(self.enc_concentration_levels_assertion(self.current_level+1))
+            
+            print("\n{:-^70}".format("[ Working at level=" + str(self.current_level) + " ]"))
+            stdout.flush()
+
+            # reachability test:
+            print("[" + colour_str(C_BOLD, "i") + "] Adding the reachability test...")       
+            self.solver.push()
+
+            self.solver_add(self.enc_state_with_blocking(self.current_level,state))
+                
+            result = self.solver.check()
+            if result == sat:
+                print("[" + colour_str(C_BOLD, "+") + "] " + colour_str(C_GREEN, "SAT at level=" + str(self.current_level)))
+                if print_witness:
+                    print("\n{:=^70}".format("[ WITNESS ]"))
+                    self.decode_witness(self.current_level)
+                break
+            else:
+                self.solver.pop()
+
+            print("[" + colour_str(C_BOLD, "i") + "] Unrolling the transition relation")
+            self.solver_add(self.enc_transition_relation(self.current_level))
+
+            print("{:->70}".format("[ level=" + str(self.current_level) + " done ]"))
+            self.current_level += 1
+
+            if self.current_level > max_level:
+                print("Stopping at level=" + str(max_level))
+                break
+
+        if print_time:
+            # stop = time()
+            stop = resource.getrusage(resource.RUSAGE_SELF).ru_utime
+            self.verification_time = stop-start
+            print()
+            print("\n[i] {: >60}".format(" Time: " + repr(self.verification_time) + " s"))
+            
+        if print_mem:
+            print("[i] {: >60}".format(" Memory: " + repr(resource.getrusage(resource.RUSAGE_SELF).ru_maxrss/(1024*1024)) + " MB"))
+    
+    
+    def get_verification_time(self):
+        return self.verification_time
+    
+    
+    def show_encoding(self, state, print_witness=True, 
+            print_time=False, print_mem=False, max_level=100):
+        """Encoding debug function"""
+
+        self.reset()
+
+        self.prepare_all_variables()
+        init_s = self.enc_init_state(0)
+        print(init_s)
+        self.solver_add(init_s)
+        self.current_level = 0
+
+        self.prepare_all_variables()
+        
+        while True:
+            self.prepare_all_variables()
+
+            print("-----[ Working at level=" + str(self.current_level) + " ]-----")
+            stdout.flush()
+
+            # reachability test:
+            print("[i] Adding the reachability test...")       
+            self.solver.push()
+
+            s = self.enc_min_state(self.current_level,state)
+            print("Test: ", s)
+            
+            self.solver_add(s)
+
+            result = self.solver.check()
+            if result == sat:
+                print("\n[+] " + colour_str(C_RED, "SAT at level=" + str(self.current_level)))
+                if print_witness:
+                    self.decode_witness(self.current_level)
+                break
+            else:
+                self.solver.pop()
+
+            print("[i] Unrolling the transition relation")
+            t = self.enc_transition_relation(self.current_level)
+            print(t)
+            self.solver_add(t)
+
+            print("-----[ level=" + str(self.current_level) + " done ]")
+            self.current_level += 1
+
+            if self.current_level > max_level:
+                print("Stopping at level=" + str(max_level))
+                break
+            else:
+                x=input("Next level? ")
+                x=x.lower()
+                if not (x == "y" or x == "yes"):
+                    break
+
+
+# EOF