Moved files to 'reactics-smt'

reactics-smt/Makefile

@@ -0,0 +1,8 @@
+# Makefile
+
+clean:
+	rm -rf rs/__pycache__ smt/rs/__pycache__ logics/__pycache__ __pycache__
+	find . -name '*.pyc' -exec rm -f {} \;
+	
+cleanall: clean
+	rm -rf *.log

reactics-smt/colour.py

@@ -0,0 +1,33 @@
+C_HEADER = '\033[95m'
+C_BLUE = '\033[94m'
+C_GREEN = '\033[92m'
+C_SOMEOTHER = '\033[93m'
+C_RED = '\033[91m'
+C_ENDC = '\033[0m'
+C_BOLD = '\033[1m'
+C_UNDERLINE = '\033[4m'
+
+C_MARK_INFO = C_BOLD + "[" + C_GREEN + "*" + C_ENDC + C_BOLD + "]" + C_ENDC
+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

reactics-smt/doc/overview.py

@@ -0,0 +1,44 @@
+#!/usr/bin/env python
+
+from rs import *
+from smt import *
+from logics import *
+from rsltl_shortcuts import *
+
+def ex():
+    
+    prs = ReactionSystemWithConcentrationsParam()
+    
+    ent_with_conc = [("a", 3), ("b",2), ("c",1), ("h",1)]
+
+    for ec in ent_with_conc:
+        prs.add_bg_set_entity(ec)
+
+    lda = prs.get_param("lda")
+    
+    prs.add_reaction([("a",1)],[("h",1)],[("b", 2)])
+    prs.add_reaction(lda,[("h",1)],[("c",1)])
+    
+    ##
+
+    ca = ContextAutomatonWithConcentrations(prs)
+    ca.add_init_state("0")
+    ca.add_state("1")
+    
+    ca.add_transition("0", [("a", 3)], "1")
+    ca.add_transition("1", [], "1")
+    ca.add_transition("1", [("h", 1)], "1")
+
+    crprs = ReactionSystemWithAutomaton(prs, ca)
+    crprs.show()
+
+    pc = param_entity(lda, "a") == 0
+    f = ltl_F(bag_entity("h") == 0, "c")
+
+    checker = SmtCheckerRSCParam(crprs, optimise=True)
+    checker.check_rsltl(formulae_list=[f], param_constr=pc)
+    
+if __name__ == "__main__":
+	ex()
+
+# EOF

reactics-smt/logics/__init__.py

@@ -0,0 +1,5 @@
+from logics.rsltl import Formula_rsLTL
+from logics.bags import BagDescription
+from logics.param_constr import ParamConstraint
+from logics.rsltl_encoder import rsLTL_Encoder
+from logics.param_constr_encoder import ParamConstr_Encoder

reactics-smt/logics/bags.py

@@ -0,0 +1,89 @@
+from enum import Enum
+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
+        self.right_operand = R_oper
+        self.entity = entity
+
+        self.sanity_check()
+
+    def __repr__(self):
+        if self.f_type == BagDesc_oper.entity:
+            return self.entity
+        if self.f_type == BagDesc_oper.true:
+            return "TRUE"
+        if self.f_type == BagDesc_oper.l_and:
+            return "(" + repr(self.left_operand) + " & " + repr(self.right_operand) + ")"
+        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) + ")"
+        if self.f_type == BagDesc_oper.lt:
+            return repr(self.left_operand) + " < " + repr(self.right_operand)
+        if self.f_type == BagDesc_oper.le:
+            return repr(self.left_operand) + " <= " + repr(self.right_operand)
+        if self.f_type == BagDesc_oper.eq:
+            return repr(self.left_operand) + " == " + repr(self.right_operand)
+        if self.f_type == BagDesc_oper.ge:
+            return repr(self.left_operand) + " >= " + repr(self.right_operand)
+        if self.f_type == BagDesc_oper.gt:
+            return repr(self.left_operand) + " > " + repr(self.right_operand)
+
+    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: {:s} (type: {:s})".format(
+                            str(operand), str(type(operand))))
+
+    @classmethod
+    def f_entity(cls, entity_name):
+        return cls(BagDesc_oper.entity, entity=entity_name)
+
+    @classmethod
+    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)
+
+    def __le__(self, other):
+        return BagDescription(BagDesc_oper.le, L_oper=self, R_oper=other)
+
+    def __eq__(self, other):
+        return BagDescription(BagDesc_oper.eq, L_oper=self, R_oper=other)
+
+    def __ge__(self, other):
+        return BagDescription(BagDesc_oper.ge, L_oper=self, R_oper=other)
+
+    def __gt__(self, other):
+        return BagDescription(BagDesc_oper.gt, L_oper=self, R_oper=other)
+
+    def __and__(self, other):
+        return BagDescription(BagDesc_oper.l_and, L_oper=self, R_oper=other)
+
+    def __or__(self, other):
+        return BagDescription(BagDesc_oper.l_or, L_oper=self, R_oper=other)
+
+    def __invert__(self):
+        return BagDescription(BagDesc_oper.l_not, L_oper=self)
+
+# EOF

reactics-smt/logics/param_constr.py

@@ -0,0 +1,99 @@
+from enum import Enum
+ParamConstraint_oper = Enum(
+    'ParamConstraint_oper',
+    'param_entity true l_and l_or l_not lt le eq ge gt')
+
+
+class ParamConstraint(object):
+
+    def __init__(self, f_type, L_oper=None, R_oper=None, param=None, entity=""):
+        self.f_type = f_type
+        self.left_operand = L_oper
+        self.right_operand = R_oper
+        self.param = param
+        self.entity = entity
+
+        self.sanity_check()
+
+    def __repr__(self):
+        if self.f_type == ParamConstraint_oper.param_entity:
+            return "{:s}[{:s}]".format(self.param.name, self.entity)
+        if self.f_type == ParamConstraint_oper.true:
+            return "TRUE"
+        if self.f_type == ParamConstraint_oper.l_and:
+            return "(" + repr(self.left_operand) + " & " + repr(self.right_operand) + ")"
+        if self.f_type == ParamConstraint_oper.l_or:
+            return "(" + repr(self.left_operand) + " | " + repr(self.right_operand) + ")"
+        if self.f_type == ParamConstraint_oper.l_not:
+            return "~(" + repr(self.left_operand) + ")"
+        if self.f_type == ParamConstraint_oper.lt:
+            return repr(self.left_operand) + " < " + repr(self.right_operand)
+        if self.f_type == ParamConstraint_oper.le:
+            return repr(self.left_operand) + " <= " + repr(self.right_operand)
+        if self.f_type == ParamConstraint_oper.eq:
+            return repr(self.left_operand) + " == " + repr(self.right_operand)
+        if self.f_type == ParamConstraint_oper.ge:
+            return repr(self.left_operand) + " >= " + repr(self.right_operand)
+        if self.f_type == ParamConstraint_oper.gt:
+            return repr(self.left_operand) + " > " + repr(self.right_operand)
+
+    def sanity_check(self):
+        """Sanity checks"""
+        for operand in (self.left_operand, self.right_operand):
+            if operand:
+                if not(
+                        isinstance(operand, ParamConstraint)
+                        or isinstance(operand, int)):
+                    raise RuntimeError(
+                        "Unexpected operand type for a bag: {:s} (type: {:s})".format(
+                            str(operand), str(type(operand))))
+
+    @classmethod
+    def f_param_ent(cls, param, entity_name):
+        return cls(ParamConstraint_oper.param_entity, param=param,
+                   entity=entity_name)
+
+    @classmethod
+    def f_TRUE(cls):
+        return cls(ParamConstraint_oper.true)
+
+    @classmethod
+    def f_And(cls, arg_L, arg_R):
+        return cls(ParamConstraint_oper.l_and, L_oper=arg_L, R_oper=arg_R)
+
+    @classmethod
+    def f_Not(cls, arg):
+        return cls(ParamConstraint_oper.l_not, L_oper=arg)
+
+    def __lt__(self, other):
+        return ParamConstraint(
+            ParamConstraint_oper.lt, L_oper=self, R_oper=other)
+
+    def __le__(self, other):
+        return ParamConstraint(
+            ParamConstraint_oper.le, L_oper=self, R_oper=other)
+
+    def __eq__(self, other):
+        return ParamConstraint(
+            ParamConstraint_oper.eq, L_oper=self, R_oper=other)
+
+    def __ge__(self, other):
+        return ParamConstraint(
+            ParamConstraint_oper.ge, L_oper=self, R_oper=other)
+
+    def __gt__(self, other):
+        return ParamConstraint(
+            ParamConstraint_oper.gt, L_oper=self, R_oper=other)
+
+    def __and__(self, other):
+        return ParamConstraint(
+            ParamConstraint_oper.l_and, L_oper=self, R_oper=other)
+
+    def __or__(self, other):
+        return ParamConstraint(
+            ParamConstraint_oper.l_or, L_oper=self, R_oper=other)
+
+    def __invert__(self):
+        return ParamConstraint(ParamConstraint_oper.l_not, L_oper=self)
+
+# EOF

reactics-smt/logics/param_constr_encoder.py

@@ -0,0 +1,70 @@
+from logics.param_constr import *
+from z3 import And, Not, Or
+
+
+class ParamConstr_Encoder(object):
+    """Class for encoding parameter constraints"""
+
+    def __init__(self, smt_checker):
+        self.smt_checker = smt_checker
+        self.rs = smt_checker.rs
+
+        # self.v = None
+        # self.v_ctx = None
+        # self.loop_position = None
+
+    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 encode(self, param_constr):
+
+        if not param_constr:
+            raise RuntimeError("param_constr is None")
+
+        if param_constr.f_type == ParamConstraint_oper.param_entity:
+            return self.smt_checker.get_enc_param(
+                param_constr.param.name, param_constr.entity)
+
+        if param_constr.f_type == ParamConstraint_oper.true:
+            return True
+
+        if param_constr.f_type == ParamConstraint_oper.l_and:
+            return And(self.encode(param_constr.left_operand),
+                       self.encode(param_constr.right_operand))
+
+        if param_constr.f_type == ParamConstraint_oper.l_or:
+            return Or(self.encode(param_constr.left_operand),
+                      self.encode(param_constr.right_operand))
+
+        if param_constr.f_type == ParamConstraint_oper.l_not:
+            return Not(self.encode(param_constr.left_operand))
+
+        if param_constr.f_type == ParamConstraint_oper.lt:
+            return self.encode(
+                param_constr.left_operand) < int(
+                param_constr.right_operand)
+
+        if param_constr.f_type == ParamConstraint_oper.le:
+            return self.encode(
+                param_constr.left_operand) <= int(
+                param_constr.right_operand)
+
+        if param_constr.f_type == ParamConstraint_oper.eq:
+            return self.encode(
+                param_constr.left_operand) == int(
+                param_constr.right_operand)
+
+        if param_constr.f_type == ParamConstraint_oper.ge:
+            return self.encode(
+                param_constr.left_operand) >= int(
+                param_constr.right_operand)
+
+        if param_constr.f_type == ParamConstraint_oper.gt:
+            return self.encode(
+                param_constr.left_operand) > int(
+                param_constr.right_operand)
+
+        assert False, "Unsupported case {:s}".format(param_constr.f_type)

reactics-smt/logics/rsltl.py

@@ -0,0 +1,119 @@
+from z3 import *
+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')
+
+
+class Formula_rsLTL(object):
+
+    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)
+        if self.f_type == rsLTL_form_type.l_not:
+            return "~(" + repr(self.left_operand) + ")"
+        if self.f_type == rsLTL_form_type.t_globally:
+            return "G[" + repr(self.sub_operand) + "](" + repr(self.left_operand) + ")"
+        if self.f_type == rsLTL_form_type.t_finally:
+            return "F[" + repr(self.sub_operand) + "](" + repr(self.left_operand) + ")"
+        if self.f_type == rsLTL_form_type.t_next:
+            return "X[" + repr(self.sub_operand) + "](" + repr(self.left_operand) + ")"
+        if self.f_type == rsLTL_form_type.l_and:
+            return "(" + repr(self.left_operand) + " & " + repr(self.right_operand) + ")"
+        if self.f_type == rsLTL_form_type.l_or:
+            return "(" + repr(self.left_operand) + " | " + repr(self.right_operand) + ")"
+        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) + ")"
+        if self.f_type == rsLTL_form_type.t_release:
+            return "(" + repr(
+                self.left_operand) + " R[" + repr(
+                self.sub_operand) + "] " + repr(
+                self.right_operand) + ")"
+
+    @property
+    def is_bag(self):
+        return self.f_type == rsLTL_form_type.bag
+
+    @classmethod
+    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)
+
+    @classmethod
+    def f_Or(cls, arg_L, arg_R):
+        return cls(rsLTL_form_type.l_or, L_oper=arg_L, R_oper=arg_R)
+
+    @classmethod
+    def f_Implies(cls, arg_L, arg_R):
+        return cls(rsLTL_form_type.l_implies, L_oper=arg_L, R_oper=arg_R)
+
+    @classmethod
+    def f_X(cls, sub, arg):
+        return cls(rsLTL_form_type.t_next, L_oper=arg, sub_oper=sub)
+
+    @classmethod
+    def f_G(cls, sub, arg):
+        return cls(rsLTL_form_type.t_globally, L_oper=arg, sub_oper=sub)
+
+    @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)
+
+    @classmethod
+    def f_F(cls, sub, arg_L):
+        return cls(rsLTL_form_type.t_finally, L_oper=arg_L, sub_oper=sub)
+
+    @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)
+
+    def __and__(self, other):
+        return Formula_rsLTL(rsLTL_form_type.l_and, L_oper=self, R_oper=other)
+
+    def __or__(self, other):
+        return Formula_rsLTL(rsLTL_form_type.l_or, L_oper=self, R_oper=other)
+
+    def __invert__(self):
+        return Formula_rsLTL(rsLTL_form_type.l_not, L_oper=self)
+
+
+# EOF

reactics-smt/logics/rsltl_encoder.py

@@ -0,0 +1,400 @@
+from logics.rsltl import *
+
+
+def simplify(x):
+    return x
+
+
+class rsLTL_Encoder(object):
+    """Class for encoding rsLTL formulae for a given smt_checker instance"""
+
+    def __init__(self, smt_checker):
+        self.smt_checker = smt_checker
+        self.rs = smt_checker.rs
+
+        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):
+        self.ncalls_encode = 0
+        self.ncalls_encode_approx = 0
+
+    def cache_init(self, bound):
+        """
+        Cache for formulae encodings
+        """
+        self.cache_hits = 0
+        self.enc_fcache = [{} for level in range(0, bound+1)]
+        self.enc_fcache_approx = [{} for level in range(0, bound+1)]
+
+    def cache_save(self, formula, level, formula_encoding):
+        self.enc_fcache[level][formula] = formula_encoding
+
+    def cache_approx_save(self, formula, level, formula_encoding):
+        self.enc_fcache_approx[level][formula] = formula_encoding
+
+    def cache_query(self, formula, level):
+        if formula in self.enc_fcache[level]:
+            self.cache_hits += 1
+            return self.enc_fcache[level][formula]
+        else:
+            return None
+
+    def cache_query_approx(self, formula, level):
+        if formula in self.enc_fcache_approx[level]:
+            self.cache_hits += 1
+            return self.enc_fcache_approx[level][formula]
+        else:
+            return None
+
+    def get_cache_hits(self):
+        return self.cache_hits
+
+    def flush_cache(self):
+        self.enc_fcache = None
+        self.enc_fcache_approx = None
+
+    def get_ncalls(self):
+        return (self.ncalls_encode, self.ncalls_encode_approx)
+
+    def encode_bag(self, bag_formula, level, context=False):
+
+        if not bag_formula:
+            raise RuntimeError("bag_formula is None")
+
+        if bag_formula.f_type == BagDesc_oper.entity:
+            entity_id = self.rs.get_entity_id(bag_formula.entity)
+            if context:
+                return self.v_ctx[level][entity_id]
+            else:
+                return self.v[level][entity_id]
+
+        if bag_formula.f_type == BagDesc_oper.true:
+            return True
+
+        if bag_formula.f_type == BagDesc_oper.l_and:
+            return And(
+                self.encode_bag(bag_formula.left_operand, level, context),
+                self.encode_bag(bag_formula.right_operand, level, context))
+
+        if bag_formula.f_type == BagDesc_oper.l_or:
+            return Or(
+                self.encode_bag(bag_formula.left_operand, level, context),
+                self.encode_bag(bag_formula.right_operand, level, context))
+
+        if bag_formula.f_type == BagDesc_oper.l_not:
+            return Not(self.encode_bag(
+                bag_formula.left_operand, level, context))
+
+        if bag_formula.f_type == BagDesc_oper.lt:
+            return self.encode_bag(
+                bag_formula.left_operand, level, context) < int(
+                bag_formula.right_operand)
+
+        if bag_formula.f_type == BagDesc_oper.le:
+            return self.encode_bag(
+                bag_formula.left_operand, level, context) <= int(
+                bag_formula.right_operand)
+
+        if bag_formula.f_type == BagDesc_oper.eq:
+            return self.encode_bag(
+                bag_formula.left_operand, level, context) == int(
+                bag_formula.right_operand)
+
+        if bag_formula.f_type == BagDesc_oper.ge:
+            return self.encode_bag(
+                bag_formula.left_operand, level, context) >= int(
+                bag_formula.right_operand)
+
+        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):
+        res = self.encode_bag(bag_formula, level)
+        assert res is not None
+        return res
+
+    def encode_bag_ctx(self, bag_formula, level):
+        res = self.encode_bag(bag_formula, level, context=True)
+        assert res is not None
+        return res
+
+    def encode(self, formula, level, bound):
+
+        self.ncalls_encode += 1
+
+        from_cache = self.cache_query(formula, level)
+        if from_cache is not None:
+            return from_cache
+
+        enc = None
+
+        if not isinstance(formula, Formula_rsLTL):
+            raise NotImplementedError(
+                "Unsupported formula type: " + str(type(formula)))
+
+        if level > bound:
+            raise RuntimeError(
+                "level > bound. Unexpected behaviour. The encoding does not support levels higher than a bound.")
+
+        if formula.f_type == rsLTL_form_type.bag:
+            enc = self.encode_bag_state(formula.bag_descr, level)
+
+        elif formula.f_type == rsLTL_form_type.l_not:
+            subform = formula.left_operand
+            if subform.is_bag:
+                enc = Not(self.encode_bag_state(subform, level))
+            else:
+                raise RuntimeError("Negation can be applied to bags only")
+
+        elif formula.f_type == rsLTL_form_type.l_and:
+            enc = And(
+                self.encode(formula.left_operand, level, bound),
+                self.encode(formula.right_operand, level, bound)
+            )
+
+        elif formula.f_type == rsLTL_form_type.l_or:
+            enc = Or(
+                self.encode(formula.left_operand, level, bound),
+                self.encode(formula.right_operand, level, bound)
+            )
+
+        elif formula.f_type == rsLTL_form_type.l_implies:
+            enc = Implies(
+                self.encode(formula.left_operand, level, bound),
+                self.encode(formula.right_operand, level, bound)
+            )
+
+        elif formula.f_type == rsLTL_form_type.t_next:
+            if level < bound:
+                enc = And(
+                    self.encode(formula.left_operand, level + 1, bound),
+                    self.encode_bag_ctx(formula.sub_operand, level)
+                )
+            else:
+                # level == bound
+                enc = False
+                for loop_level in range(1, bound+1):
+                    enc = simplify(Or(enc, And(self.loop_position == loop_level, self.encode(
+                        formula.left_operand, loop_level, bound))))
+                enc = And(enc, self.encode_bag_ctx(formula.sub_operand, level))
+                enc = simplify(enc)
+
+        elif formula.f_type == rsLTL_form_type.t_globally:
+            if level < bound:
+                enc = And(
+                    self.encode(formula.left_operand, level, bound),
+                    self.encode_bag_ctx(formula.sub_operand, level),
+                    self.encode(formula, level + 1, bound)
+                )
+            else:
+                # level == bound
+                enc_loops = False
+                for loop_level in range(1, bound+1):
+                    enc_loops = simplify(Or(enc_loops,
+                                            And(
+                                                self.loop_position == loop_level,
+                                                self.encode_approx(formula, loop_level, bound),
+                                            )
+                                            ))
+                enc = And(
+                    self.encode(formula.left_operand, bound, bound),
+                    enc_loops,
+                    self.encode_bag_ctx(formula.sub_operand, level)
+                )
+                enc = simplify(enc)
+
+        elif formula.f_type == rsLTL_form_type.t_finally:
+            if level < bound:
+                enc = Or(
+                    self.encode(formula.left_operand, level, bound),
+                    And(
+                        self.encode_bag_ctx(formula.sub_operand, level),
+                        self.encode(formula, level + 1, bound)
+                    )
+                )
+            else:
+                # level == bound
+                enc_loops = False
+                for loop_level in range(1, bound+1):
+                    enc_loops = simplify(Or(enc_loops,
+                                            And(
+                                                self.loop_position == loop_level,
+                                                self.encode_approx(formula, loop_level, bound),
+                                            )
+                                            ))
+                    # print(enc)
+                enc = Or(self.encode(formula.left_operand, bound, bound),
+                         And(
+                    enc_loops,
+                    self.encode_bag_ctx(formula.sub_operand, level)
+                )
+                )
+
+                enc = simplify(enc)
+
+        elif formula.f_type == rsLTL_form_type.t_until:
+            if level < bound:
+                inner_enc = self.encode(formula, level + 1, bound)
+            else:
+                # level == bound
+                inner_enc = False
+                for loop_level in range(1, bound+1):
+                    inner_enc = simplify(Or(inner_enc,
+                                            And(
+                                                self.loop_position == loop_level,
+                                                self.encode_approx(formula, loop_level, bound)
+                                            )
+                                            ))
+
+            enc = Or(
+                self.encode(formula.right_operand, level, bound),
+                And(
+                    self.encode(formula.left_operand, level, bound),
+                    inner_enc,
+                    self.encode_bag_ctx(formula.sub_operand, level)
+                )
+            )
+
+        elif formula.f_type == rsLTL_form_type.t_release:
+            if level < bound:
+                inner_enc = self.encode(formula, level + 1, bound)
+            else:
+                # level == bound
+                inner_enc = False
+                for loop_level in range(1, bound+1):
+                    inner_enc = simplify(Or(inner_enc,
+                                            And(
+                                                self.loop_position == loop_level,
+                                                self.encode_approx(formula, loop_level, bound)
+                                            )
+                                            ))
+
+            enc = And(
+                self.encode(formula.right_operand, level, bound),
+                Or(
+                    self.encode(formula.left_operand, level, bound),
+                    And(
+                        inner_enc,
+                        self.encode_bag_ctx(formula.sub_operand, level)
+                    )
+                )
+            )
+
+        else:
+            raise NotImplementedError("Unsupported operator")
+
+        if enc is None:
+            raise RuntimeError("Encoding is NONE. Should never happen")
+
+        self.cache_save(formula, level, enc)
+
+        return enc
+
+    def encode_approx(self, formula, level, bound):
+        """Provides the approximation-encoding
+
+        Used by encode()
+        """
+
+        self.ncalls_encode_approx += 1
+
+        enc = None
+
+        from_cache = self.cache_query_approx(formula, level)
+        if from_cache is not None:
+            return from_cache
+
+        if formula.f_type == rsLTL_form_type.t_until:
+            if level < bound:
+                enc = Or(
+                    self.encode(formula.right_operand, level, bound),
+                    And(
+                        self.encode(formula.left_operand, level, bound),
+                        self.encode_approx(formula, level + 1, bound),
+                        self.encode_bag_ctx(formula.sub_operand, level)
+                    )
+                )
+            else:
+                # level == bound
+                enc = self.encode(formula.right_operand, bound, bound)
+
+        elif formula.f_type == rsLTL_form_type.t_release:
+            if level < bound:
+                enc = And(
+                    self.encode(formula.right_operand, level, bound),
+                    Or(
+                        self.encode(formula.left_operand, level, bound),
+                        And(
+                            self.encode_approx(formula, level + 1, bound),
+                            self.encode_bag_ctx(formula.sub_operand, level)
+                        )
+                    )
+                )
+            else:
+                # level == bound
+                enc = self.encode(formula.right_operand, bound, bound)
+
+        elif formula.f_type == rsLTL_form_type.t_globally:
+            if level < bound:
+                enc = And(
+                    self.encode(formula.left_operand, level, bound),
+                    self.encode_bag_ctx(formula.sub_operand, level),
+                    self.encode_approx(formula, level + 1, bound)
+                )
+            else:
+                # level == bound
+                enc = self.encode(formula.left_operand, bound, bound)
+
+        elif formula.f_type == rsLTL_form_type.t_finally:
+            if level < bound:
+                enc = Or(
+                    self.encode(formula.left_operand, level, bound),
+                    And(
+                        self.encode_bag_ctx(formula.sub_operand, level),
+                        self.encode_approx(formula, level + 1, bound)
+                    )
+                )
+            else:
+                # level == bound
+                enc = self.encode(formula.left_operand, bound, bound)
+
+        else:
+            raise NotImplementedError(
+                "Unsupported operator in approximation encoding")
+
+        if enc is None:
+            raise RuntimeError("Encoding is NONE. Should never happen")
+
+        self.cache_approx_save(formula, level, enc)
+
+        return enc

reactics-smt/rs/__init__.py

@@ -0,0 +1,14 @@
+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
+
+from rs.network_of_context_automata import NetworkOfContextAutomata
+from rs.reaction_system_with_automaton import ReactionSystemWithAutomaton
+from rs.reaction_system_with_autnet import ReactionSystemWithNetworkOfAutomata
+
+# EOF

reactics-smt/rs/context_automaton.py

@@ -0,0 +1,171 @@
+from sys import exit
+from colour import *
+
+
+class ContextAutomaton(object):
+
+    def __init__(self, reaction_system):
+        self._states = []
+        self._transitions = []
+        self._init_state = None
+        self._reaction_system = reaction_system
+        self._name = ""
+        self._prod_entities = set()
+
+    @property
+    def states(self):
+        return self._states
+
+    @property
+    def transitions(self):
+        return self._transitions
+
+    @property
+    def prod_entities(self):
+        return self._prod_entities
+
+    @property
+    def reaction_system(self):
+        return self._reaction_system
+
+    @property
+    def name(self):
+        return self._name
+
+    @name.setter
+    def name(self, automaton_name):
+        self._name = automaton_name
+
+    def add_state(self, name):
+        if name not in self._states:
+            self._states.append(name)
+        else:
+            print("\'%s\' already added. skipping..." % (name,))
+
+    def add_states(self, states_set):
+        for st in states_set:
+            self.add_state(st)
+
+    def add_init_state(self, name):
+        self.add_state(name)
+        self._init_state = self._states.index(name)
+
+    def get_init_state_name(self):
+        if self._init_state == None:
+            return None
+        return self._states[self._init_state]
+
+    def is_state(self, name):
+        if name in self._states:
+            return True
+        else:
+            return False
+
+    def get_state_id(self, name):
+        try:
+            return self._states.index(name)
+        except ValueError:
+            print_error("Undefined context automaton state: " + repr(name))
+            exit(1)
+
+    def get_state_name(self, state_id):
+        return self._states[state_id]
+
+    def get_init_state_id(self):
+        return self._init_state
+
+    def print_states(self):
+        for state in self._states:
+            print(state)
+
+    def is_valid_rs_set(self, elements):
+        if set(elements).issubset(self._reaction_system.background_set):
+            return True
+        else:
+            return False
+
+    def is_valid_context(self, context):
+        return self.is_valid_rs_set(context)
+
+    def get_set_of_ids(self, elements):
+        """Converts a set/list/tuple of entities into a set of their ids"""
+
+        new_set = set()
+        for e in set(elements):
+            new_set.add(self._reaction_system.get_entity_id(e))
+        return new_set
+
+    def add_transition(self, src, context_set, dst):
+        if not type(context_set) is set and not type(context_set) is list:
+            print("Contexts set must be of type set or list")
+
+        if not self.is_valid_context(context_set):
+            raise RuntimeError(
+                "one of the entities in the context set is unknown (undefined)!")
+
+        if not self.is_state(src):
+            raise RuntimeError(
+                "\"" + src + "\" is an unknown (undefined) state")
+
+        if not self.is_state(dst):
+            raise RuntimeError(
+                "\"" + dst + "\" is an unknown (undefined) state")
+
+        new_context_set = set()
+        for e in set(context_set):
+            new_context_set.add(self._reaction_system.get_entity_id(e))
+
+        self._prod_entities |= new_context_set
+
+        self._transitions.append(
+            (self.get_state_id(src),
+             new_context_set, self.get_state_id(dst)))
+
+    def rsset2str(self, elements):
+        """Converts the set of entities ids into the string with their names"""
+        if len(elements) == 0:
+            return "0"
+        s = "{"
+        for c in elements:
+            s += " " + self._reaction_system.get_entity_name(c)
+        s += " }"
+        return s
+
+    def context2str(self, ctx):
+        return self.rsset2str(ctx)
+
+    def show_transitions(self):
+        print(C_MARK_INFO + " Context automaton transitions:")
+        for transition in self._transitions:
+            str_transition = str(transition[0]) + " --( "
+            str_transition += self.context2str(transition[1])
+            str_transition += " )--> " + str(transition[2])
+            print(" - " + str_transition)
+
+    def show_states(self):
+        init_state_name = self.get_init_state_name()
+        print(C_MARK_INFO + " Context automaton states:")
+        for state in self._states:
+            print(" - " + state, end="")
+            if state == init_state_name:
+                print(" [init]")
+            else:
+                print()
+
+    def show_header(self):
+        if self.name:
+            name_string = ": " + colour_str(C_BOLD, self.name)
+            print(C_MARK_INFO + " Context automaton" + name_string)
+
+    def show_prod_entities(self):
+        print(C_MARK_INFO + " Context automaton possible products:")
+        for entity in self._prod_entities:
+            print(" - " + self._reaction_system.get_entity_name(entity))
+
+    def show(self):
+        self.show_header()
+        self.show_states()
+        self.show_transitions()
+        self.show_prod_entities()
+
+# EOF

reactics-smt/rs/context_automaton_with_concentrations.py

@@ -0,0 +1,83 @@
+from sys import exit
+from colour import *
+
+from rs.context_automaton import ContextAutomaton
+
+
+class ContextAutomatonWithConcentrations(ContextAutomaton):
+
+    def __init__(self, reaction_system):
+        super(ContextAutomatonWithConcentrations,
+              self).__init__(reaction_system)
+
+    def is_valid_context(self, context):
+        if set(
+                [e for e, lvl in context]).issubset(
+                self._reaction_system.background_set):
+            return True
+        else:
+            return False
+
+    def context2str(self, ctx):
+        if len(ctx) == 0:
+            return "0"
+        s = "{"
+        for ent, lvl in ctx:
+            s += " " + str((self._reaction_system.get_entity_name(ent), lvl))
+        s += " }"
+        return s
+
+    def add_transition(self, src, context_set, dst):
+        if not type(context_set) is set and not type(context_set) is list:
+            print("Contexts set must be of type set or list")
+
+        if not self.is_valid_context(context_set):
+            raise RuntimeError(
+                "one of the entities in the context set is unknown (undefined)!")
+
+        if not self.is_state(src):
+            raise RuntimeError(
+                "\"" + src + "\" is an unknown (undefined) state")
+
+        if not self.is_state(dst):
+            raise RuntimeError(
+                "\"" + dst + "\" is an unknown (undefined) state")
+
+        new_context_set = set()
+        for ent, lvl in set(context_set):
+            new_context_set.add(
+                (self._reaction_system.get_entity_id(ent), lvl))
+
+        self._transitions.append(
+            (self.get_state_id(src),
+             new_context_set, self.get_state_id(dst)))
+
+    def get_automaton_with_flat_contexts(self, ordinary_reaction_system):
+
+        ca = ContextAutomaton(ordinary_reaction_system)
+        ca._states = self._states
+        ca._init_state = self._init_state
+
+        for src, ctx, dst in self._transitions:
+
+            new_ctx = set()
+
+            for ent, conc in ctx:
+                for i in range(1, conc+1):
+                    n = self._reaction_system.get_entity_name(
+                        ent) + "#" + str(i)
+                    ca._reaction_system.ensure_bg_set_entity(n)
+                    new_ctx.add(n)
+
+            ca.add_transition(
+                ca.get_state_name(src),
+                new_ctx, ca.get_state_name(dst))
+
+        return ca
+
+    def show(self):
+        self.show_header()
+        self.show_states()
+        self.show_transitions()
+
+# EOF

reactics-smt/rs/extended_context_automaton.py

@@ -0,0 +1,176 @@
+from sys import exit
+from colour import *
+
+from rs.context_automaton import ContextAutomaton
+
+
+class ExtendedContextAutomaton(ContextAutomaton):
+    """Extended Context Automaton
+
+    Supports transitions with actions.
+
+    Each transitions is additionally guarded with
+    reactants and inhibitors.
+
+    The provided context entities are the products
+    of the reactions labelling the transition taken.
+    """
+
+    def __init__(self, reaction_system):
+        super(ExtendedContextAutomaton, self).__init__(reaction_system)
+        self._actions = []
+        self._transitions_for_products = dict()
+        self._actions_for_products = dict()
+
+    @property
+    def number_of_actions(self):
+        return len(self._actions)
+
+    @property
+    def actions(self):
+        return self._actions
+
+    def has_action(self, action):
+        """Checks if the automaton supports a given action"""
+
+        return action in self._actions
+
+    def get_transitions_producing_entity(self, entity):
+        """Returns the transitions that produce a given entity"""
+
+        if entity in self._transitions_for_products:
+            return self._transitions_for_products[entity]
+        else:
+            return []
+
+    def get_actions_producing_entity(self, entity):
+        """Returns the actions that produce a given entity"""
+
+        if entity in self._actions_for_products:
+            return self._actions_for_products[entity]
+        else:
+            return set()
+
+    def can_produce_entity(self, entity):
+        """Check if the automaton can produce an entity"""
+
+        if entity in self._actions_for_products:
+            return True
+        else:
+            return False
+
+    def add_transition(self, src, actions, ctx_reaction, dst):
+        """Adds a transition
+
+        src: is the source state name
+        dst: is the destination state name
+        actions: is the set of actions with which the transitions is synchronised
+        ctx_reaction: is the context reaction associated with the transition
+        """
+
+        ctx_reactants, ctx_inhibitors, ctx_products = ctx_reaction
+
+        if not type(ctx_products) is set and not type(ctx_products) is list:
+            print("Contexts set (context products) must be of type set or list")
+
+        if not self.is_valid_rs_set(ctx_reactants):
+            raise RuntimeError(
+                "one of the entities in the reactants set is unknown (undefined)!")
+
+        if not self.is_valid_rs_set(ctx_inhibitors):
+            raise RuntimeError(
+                "one of the entities in the inhibitors set is unknown (undefined)!")
+
+        if not self.is_valid_rs_set(ctx_products):
+            raise RuntimeError(
+                "one of the entities in the context set is unknown (undefined)!")
+
+        if not self.is_state(src):
+            raise RuntimeError(
+                "\"" + src + "\" is an unknown (undefined) state")
+
+        if not self.is_state(dst):
+            raise RuntimeError(
+                "\"" + dst + "\" is an unknown (undefined) state")
+
+        src_id = self.get_state_id(src)
+        dst_id = self.get_state_id(dst)
+        act_ids = self.get_set_of_action_ids(actions)
+        r_ids = self.get_set_of_ids(ctx_reactants)
+        i_ids = self.get_set_of_ids(ctx_inhibitors)
+        p_ids = self.get_set_of_ids(ctx_products)
+
+        new_transition = (src_id, act_ids, (r_ids, i_ids, p_ids), dst_id)
+
+        for product_id in p_ids:
+            self._transitions_for_products.setdefault(product_id, [])
+            self._transitions_for_products[product_id].append(new_transition)
+            self._actions_for_products.setdefault(product_id, set())
+            self._actions_for_products[product_id] |= set(actions)
+
+        self._prod_entities |= p_ids
+
+        self._transitions.append(new_transition)
+
+    def show_transitions(self):
+        """Prints the set of registered transitions"""
+
+        print(C_MARK_INFO + " Context automaton transitions:")
+        for src_id, act_id, reaction, dst_id in self._transitions:
+            str_transition = self.get_state_name(src_id) + " --( "
+            str_transition += "<" + self.get_actions_str(act_id) + "> | "
+            str_transition += "( " + self.rsset2str(reaction[0]) + "," + self.rsset2str(
+                reaction[1]) + "," + self.rsset2str(reaction[2]) + " )"
+            str_transition += " )--> " + self.get_state_name(dst_id)
+            print(" - " + str_transition)
+
+    def add_action(self, action_name):
+        """Registers an action"""
+
+        if action_name not in self._actions:
+            self._actions.append(action_name)
+        else:
+            print("\'%s\' already added. skipping..." % (action_name,))
+
+    def get_action_id(self, action_name):
+        """For an action name returns its id"""
+
+        try:
+            return self._actions.index(action_name)
+        except ValueError:
+            print_error("Undefined context automaton action: " +
+                        repr(action_name))
+            exit(1)
+
+    def get_action_name(self, action_id):
+        return self._actions[action_id]
+
+    def get_set_of_action_ids(self, actions):
+        """Converts a set of actions into the set of their ids"""
+
+        act_ids = set()
+        for act in actions:
+            act_ids.add(self.get_action_id(act))
+        return act_ids
+
+    def get_actions_str(self, actions):
+        """Returns the string for the set of action ids given by actions"""
+
+        s = ""
+        for act in actions:
+            s += self.get_action_name(act) + ", "
+        s = s[:-2]
+        return s
+
+    def show_actions(self):
+        """Prints all the actions"""
+
+        print(C_MARK_INFO + " Context automaton actions:")
+        for act in self._actions:
+            print(" - " + act + " (id=" + str(self.get_action_id(act)) + ")")
+
+    def show(self):
+        super(ExtendedContextAutomaton, self).show()
+        self.show_actions()
+
+# EOF

reactics-smt/rs/network_of_context_automata.py

@@ -0,0 +1,111 @@
+from sys import exit
+from colour import *
+
+
+class NetworkOfContextAutomata(object):
+
+    def __init__(self, reaction_system, context_automata):
+        self.automata = []
+        self._reaction_system = reaction_system
+        self._actions = set()
+        self._prod_entities = set()
+        self._automata_for_actions = dict()
+        self._actions_for_products = dict()
+
+        if len(context_automata) < 1:
+            print("Context automata network must contain at least one automaton!")
+            exit(1)
+
+        for automaton in context_automata:
+            self.add(automaton)
+
+        self.sanity_check()
+
+    @property
+    def number_of_automata(self):
+        return len(self.automata)
+
+    @property
+    def reaction_system(self):
+        return self._reaction_system
+
+    @property
+    def prod_entities(self):
+        """Returns the set of entities that can potentially be produced by the automata in the network"""
+        return self._prod_entities
+
+    @property
+    def automata_ids(self):
+        return set(range(len(self.automata)))
+
+    def sanity_check(self):
+        """Performs a sanity check of the network of automata"""
+
+        for automaton in self.automata:
+            if automaton.reaction_system != self._reaction_system:
+                print_error(
+                    "Mismatching reaction system used in \"" +
+                    str(automaton.name) + "\"!!!")
+                exit(1)
+
+    def show_prod_entities(self):
+        print(
+            C_MARK_INFO +
+            " Possible context-products for the network of automata:")
+        for entity in self.prod_entities:
+            print(" - " + self._reaction_system.get_entity_name(entity))
+
+    def show_actions(self):
+        print(C_MARK_INFO + " Actions of the network:")
+        for action in self._actions:
+            print(" - " + action)
+
+    def register_action(self, action, aut):
+        """Associates an action with an automaton"""
+
+        self._automata_for_actions.setdefault(action, set())
+        aut_index = self.automata.index(aut)
+        self._automata_for_actions[action].add(aut_index)
+
+    def get_actions_producing_entity(self, entity):
+        """Returns the set of actions producing an entity"""
+
+        if entity in self._actions_for_products:
+            return self._actions_for_products[entity]
+        else:
+            return set()
+
+    def get_automata_with_action(self, action):
+        """Returns the set of automata that support an action"""
+
+        if action in self._automata_for_actions:
+            return self._automata_for_actions[action]
+        else:
+            return set()
+
+    def add(self, aut):
+        """Adds an automaton to the network"""
+
+        self.automata.append(aut)
+        self._prod_entities |= aut.prod_entities
+        self._actions |= set(aut.actions)
+
+        for action in aut.actions:
+            self.register_action(action, aut)
+
+        for entity in aut.prod_entities:
+            self._actions_for_products.setdefault(entity, set())
+            self._actions_for_products[entity] |= aut.get_actions_producing_entity(
+                entity)
+
+    def show(self):
+        print()
+        print(C_MARK_INFO + " NETWORK OF CONTEXT AUTOMATA")
+        for ca in self.automata:
+            print()
+            ca.show()
+        print()
+        self.show_prod_entities()
+        self.show_actions()
+
+# EOF

reactics-smt/rs/reaction_system.py

@@ -0,0 +1,203 @@
+from sys import exit
+from colour import *
+
+
+class ReactionSystem(object):
+
+    def __init__(self):
+
+        self.reactions = []
+        self.background_set = []
+
+        # self.reactions_by_agents = [] # each element is 'reactions_by_prod'
+        self.reactions_by_prod = None
+
+        # legacy:
+        self.init_contexts = []
+        self.context_entities = []
+
+    @property
+    def background_set_size(self):
+        return len(self.background_set)
+
+    @property
+    def set_of_bgset_ids(self):
+        return set(range(self.background_set_size))
+
+    @property
+    def ordered_list_of_bgset_ids(self):
+        return list(range(self.background_set_size))
+
+    def assume_not_in_bgset(self, name):
+        if self.is_in_background_set(name):
+            raise RuntimeError(
+                "The entity " + name + " is already on the list")
+
+    def add_bg_set_entity(self, name):
+        self.assume_not_in_bgset(name)
+        self.background_set.append(name)
+
+    def ensure_bg_set_entity(self, name):
+        if not self.is_in_background_set(name):
+            self.background_set.append(name)
+
+    def add_bg_set_entities(self, elements):
+        for e in elements:
+            self.add_bg_set_entity(e)
+
+    def is_in_background_set(self, entity):
+        """Checks if the given name is valid wrt the background set="""
+        if entity in self.background_set:
+            return True
+        else:
+            return False
+
+    def get_entity_id(self, name):
+        try:
+            return self.background_set.index(name)
+        except ValueError:
+            print("Undefined background set entity: " + repr(name))
+            exit(1)
+
+    def get_state_ids(self, state):
+        ids = []
+        for entity in state:
+            ids.append(self.get_entity_id(entity))
+
+        return ids
+
+    def get_entity_name(self, entity_id):
+        """Returns the string corresponding to the entity"""
+        return self.background_set[entity_id]
+
+    def add_reaction(self, R, I, P):
+        """Adds a reaction"""
+
+        if R == [] or P == []:
+            raise RuntimeError("No reactants or products defined")
+
+        reactants = []
+        for entity in R:
+            reactants.append(self.get_entity_id(entity))
+
+        inhibitors = []
+        for entity in I:
+            inhibitors.append(self.get_entity_id(entity))
+
+        products = []
+        for entity in P:
+            products.append(self.get_entity_id(entity))
+
+        self.reactions.append((reactants, inhibitors, products))
+
+    def add_initial_context_set(self, context_set):
+        if context_set == []:
+            print("Empty context set is not allowed")
+            raise
+
+        integers = []
+        for entity in context_set:
+            if not entity in self.background_set:
+                print("The entity", entity, "is not in the background set")
+                raise
+            else:
+                integers.append(self.get_entity_id(entity))
+
+        self.init_contexts.append(integers)
+
+    def set_context_entities(self, entities):
+
+        for entity in entities:
+            entity_id = self.get_entity_id(entity)
+            self.context_entities.append(entity_id)
+
+    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):
+        return self.entities_ids_set_to_str(state)
+
+    def show_reactions(self, soft=False):
+        print(C_MARK_INFO + " Reactions:")
+        if soft and len(self.reactions) > 50:
+            print(" -> there are more than 50 reactions (" +
+                  str(len(self.reactions)) + ")")
+        else:
+            print(
+                " "*4 + "{0: ^35}{1: ^25}{2: ^15}".format("reactants", " inhibitors", " products"))
+            for reaction in self.reactions:
+                # print("\t( R={" + self.state_to_str(reaction[0]) + "}, I={" + self.state_to_str(reaction[1]) + "}, P={" + self.state_to_str(reaction[2]) + "} )")
+                print(" " + "- {0: ^35}{1: ^25}{2: ^15}".format("{ " + self.state_to_str(reaction[0]) + " }",
+                                                                " { " + self.state_to_str(reaction[1]) + " }",
+                                                                " { " + self.state_to_str(reaction[2]) + " }"))
+
+    def show_background_set(self):
+        print(
+            C_MARK_INFO + " Background set: {" + self.entities_names_set_to_str(self.background_set) + "}")
+
+    def show_initial_contexts(self):
+        if len(self.init_contexts) > 0:
+            print(C_MARK_INFO + " Initial context sets:")
+            for ctx in self.init_contexts:
+                print(" - {" + self.entities_ids_set_to_str(ctx) + "}")
+
+    def show_context_entities(self):
+        if len(self.context_entities) > 0:
+            print(
+                C_MARK_INFO + " Context entities: " + self.entities_ids_set_to_str(self.context_entities))
+
+    def show(self, soft=False):
+
+        self.show_background_set()
+        self.show_initial_contexts()
+        self.show_reactions(soft)
+        self.show_context_entities()
+
+    def get_reactions_by_product(self):
+        """Sorts reactions by their products and returns a dictionary of products"""
+
+        if self.reactions_by_prod != None:
+            return self.reactions_by_prod
+
+        producible_entities = set()
+
+        for reaction in self.reactions:
+            producible_entities = producible_entities.union(set(reaction[2]))
+
+        reactions_by_prod = {}
+
+        for prod_entity in producible_entities:
+            reactions_by_prod[prod_entity] = []
+            for reaction in self.reactions:
+                if prod_entity in reaction[2]:
+                    reactions_by_prod[prod_entity].append(
+                        [reaction[0], reaction[1]])
+
+        # save in cache
+        self.reactions_by_prod = reactions_by_prod
+
+        return reactions_by_prod
+
+    def sanity_check(self):
+        """Performs a sanity check on the defined reaction system"""
+
+        if self.reactions == []:
+            print("No reactions defined")
+            exit(1)
+
+        if self.background_set == []:
+            print("Empty background set")
+            exit(1)
+
+# EOF

reactics-smt/rs/reaction_system_with_autnet.py

@@ -0,0 +1,12 @@
+
+class ReactionSystemWithNetworkOfAutomata(object):
+
+    def __init__(self, reaction_system, context_automata):
+        self.rs = reaction_system
+        self.cas = context_automata
+
+    def show(self, soft=False):
+        self.rs.show(soft)
+        self.cas.show()
+
+# EOF

reactics-smt/rs/reaction_system_with_automaton.py

@@ -0,0 +1,45 @@
+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):
+
+    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):
+        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)

reactics-smt/rs/reaction_system_with_concentrations.py

@@ -0,0 +1,433 @@
+from sys import exit
+from colour import *
+
+from rs.reaction_system import ReactionSystem
+
+
+class ReactionSystemWithConcentrations(ReactionSystem):
+
+    def __init__(self):
+
+        self.reactions = []
+        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_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]
+
+    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 = []
+        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 = []
+        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 = []
+        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 add_reaction(self, R, I, P):
+        """Adds a reaction"""
+
+        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):
+        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 optimized translation to RS):")
+        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_permanent_entities()
+        self.show_meta_reactions()
+        self.show_max_concentrations()
+
+    def get_reactions_by_product(self):
+        """Sorts reactions by their products and returns a dictionary of products"""
+
+        if self.reactions_by_prod != None:
+            return self.reactions_by_prod
+
+        producible_entities = set()
+
+        for reaction in self.reactions:
+            product_entities = [e for e, c in reaction[2]]
+            producible_entities = producible_entities.union(
+                set(product_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
+                        # we append (to the end)
+                        rcts_for_p_e.append((reactants, inhibitors, products))
+                    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):
+
+        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_with_concentrations(self):
+        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)
+
+
+# EOF

reactics-smt/rs/reaction_system_with_concentrations_param.py

@@ -0,0 +1,439 @@
+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 = 1
+        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 [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_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
+
+# EOF

reactics-smt/rs_examples.py

@@ -0,0 +1,190 @@
+#!/usr/bin/env python
+
+from rs import *
+from smt import *
+import sys
+import resource
+
+def chain_reaction(print_system=False):
+    
+    if len(sys.argv) < 1+3:
+        print("provide N M B")
+        print(" B=1 - RSC")
+        print(" B=0 - Translated RSC into RS")
+        exit(1)
+
+    chainLen=int(sys.argv[1]) # chain length
+    maxConc=int(sys.argv[2]) # depth (max concentration)
+    verify_rsc=bool(int(sys.argv[3]))
+        
+    if chainLen < 1 or maxConc < 1:
+        print("be reasonable")
+        exit(1)
+    
+    r = ReactionSystemWithConcentrations()
+    r.add_bg_set_entity(("inc",1))
+    # r.add_bg_set_entity(("reset",1))
+    r.add_bg_set_entity(("dec",1))
+    # r.add_bg_set_entity(("x",1))
+    
+    for i in range(1,chainLen+1):
+        r.add_bg_set_entity(("e_" + str(i),maxConc))
+    
+    for i in range(1,chainLen+1):
+        ent = "e_" + str(i)
+        # r.add_reaction_inc(ent, "inc", [(ent, 1)],[("reset",1),(ent,maxConc)])
+        r.add_reaction_inc(ent, "inc", [(ent, 1)],[(ent,maxConc)])
+        r.add_reaction_dec(ent, "dec", [(ent, 1)],[])
+        # r.add_reaction([(ent,1),("reset",1)],[],[(ent,1)])
+        if i < chainLen:
+            r.add_reaction([(ent,maxConc)],[],[("e_"+str(i+1),1)])
+
+    r.add_reaction([("e_" + str(chainLen),maxConc)],[("dec",1)],[("e_" + str(chainLen),maxConc)])
+    # r.add_reaction([("e_" + str(chainLen),maxConc)],[],[("e_" + str(chainLen),maxConc)])
+    
+    
+    c = ContextAutomatonWithConcentrations(r)
+    c.add_init_state("init")
+    c.add_state("working")
+    c.add_transition("init", [("e_1",1),("inc",1)], "working")
+    c.add_transition("working", [("inc",1)], "working")
+    # c.add_transition("working", [("reset",1)], "working")
+    # c.add_transition("working",[("x",1)],"working")
+
+    rc = ReactionSystemWithAutomaton(r,c)
+    
+    if print_system:
+        rc.show()
+    
+    if verify_rsc:
+        smt_rsc = SmtCheckerRSC(rc)
+        prop = [('e_'+str(chainLen),maxConc)]
+        smt_rsc.check_reachability((prop,[]),max_level=maxConc*chainLen+10)
+        # smt_rsc.show_encoding(prop,print_time=True,max_level=maxConc*chainLen+10)
+
+    else:
+        orc = rc.get_ordinary_reaction_system_with_automaton()
+        if print_system:
+            print("\nTranslated:")
+            orc.show()
+        smt_tr_rs = SmtCheckerRS(orc)
+        smt_tr_rs.check_reachability(['e_'+str(chainLen)+"#"+str(maxConc)])
+    
+    # print("Reaction System with Concentrations:", smt_rsc.get_verification_time())
+    # print("Reaction System from translating RSC:", smt_tr_rs.get_verification_time())
+    
+    time=0
+    mem_usage=resource.getrusage(resource.RUSAGE_SELF).ru_maxrss/(1024*1024)
+    if verify_rsc:
+        filename_t="bench_rsc_time.log"
+        filename_m="bench_rsc_mem.log"
+        time=smt_rsc.get_verification_time()
+    else:
+        filename_t="bench_tr_rs_time.log"
+        filename_m="bench_tr_rs_mem.log"
+        time=smt_tr_rs.get_verification_time()
+
+    f=open(filename_t, 'a')    
+    log_str="(" + str(chainLen) + "," + str(maxConc) + "," + str(time) + ")\n"
+    f.write(log_str)
+    f.close()
+    
+    f=open(filename_m, 'a')
+    log_str="(" + str(chainLen) + "," + str(maxConc) + "," + str(mem_usage) + ")\n"
+    f.write(log_str)
+    f.close()
+
+
+def heat_shock_response(print_system=True,verify_rsc=True):
+
+    if len(sys.argv) < 1+1:
+        print("provide B")
+        print(" B=1 - RSC")
+        print(" B=0 - Translated RSC into RS")
+        exit(1)
+
+    verify_rsc=bool(int(sys.argv[1]))
+
+    stress_temp = 42
+    max_temp = 50
+    
+    r = ReactionSystemWithConcentrations()    
+    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(("temp",max_temp))
+    r.add_bg_set_entity(("heat",1))
+    r.add_bg_set_entity(("cool",1))
+    
+    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),("temp",stress_temp)],[],                 [("hsf",1),("hsp",1)])
+    r.add_reaction([("hsp:hsf",1)],                     [("temp",stress_temp)],[("hsp:hsf",1)])
+    r.add_reaction([("hsp",1),("hsf3:hse",1)],          [("mfp",1)],        [("hse",1),("hsp:hsf",1)])
+    r.add_reaction([("temp",stress_temp),("prot",1)],   [],                     [("mfp",1),("prot",1)])
+    r.add_reaction([("prot",1)],                        [("temp",stress_temp)], [("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)])
+    #r.add_permanency("temp",[])
+
+    c = ContextAutomatonWithConcentrations(r)
+    c.add_init_state("0")
+    c.add_state("1")
+    # c.add_transition("0", [("prot",1), ("temp",35)], "1")
+    c.add_transition("0", [("hsf",1),("prot",1),("hse",1),("temp",35)], "1")
+    #-> c.add_transition("0", [("hse",1),("prot",1),("hsp:hsf",1),("temp",stress_temp)], "1")
+    # c.add_transition("0", [("hsp",1),("prot",1),("hsf3:hse",1),("mfp",1),("hsp:mfp",1),("temp",30)], "1")
+    c.add_transition("1", [("cool",1)], "1")
+    c.add_transition("1", [("heat",1)], "1")
+    c.add_transition("1", [], "1")
+    # c.add_transition("1", [("sugar",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))
+    
+    if verify_rsc:
+        smt_rsc = SmtCheckerRSC(rc)
+        smt_rsc.check_reachability(prop,max_level=40)
+    else:
+        orc = rc.get_ordinary_reaction_system_with_automaton()
+        if print_system:
+            print("\nTranslated:")
+            orc.show()
+        smt_tr_rs = SmtCheckerRS(orc)
+        smt_tr_rs.check_reachability(rs_prop)
+
+def state_translate_rsc2rs(p):
+    return [e[0] + "#" + str(e[1]) for e in p]
+        
+# EOF

reactics-smt/rsltl_shortcuts.py

@@ -0,0 +1,108 @@
+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):
+    """
+    Assumes 0 concentration level for all the
+    not listed entities but present in 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_U(ctx_arg, a0, a1):
+    a0 = get_bag_if_str(a0)
+    a1 = get_bag_if_str(a1)
+    return Formula_rsLTL.f_U(ctx_arg, a0, a1)
+
+
+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)
+
+
+def param_entity(param, entity_name):
+    return ParamConstraint.f_param_ent(param, entity_name)
+
+
+def param_And(*args):
+    assert len(args) > 1
+    last = args[0]
+    for arg in args[1:]:
+        last = ParamConstraint.f_And(last, arg)
+    return last

reactics-smt/rssmt.py

@@ -0,0 +1,81 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+"""
+    Reaction Systems SMT-Based Model Checking Module
+
+"""
+
+import argparse
+
+from rs import *
+from smt import *
+import sys
+import rs_examples
+import rs_testing
+
+from colour import *
+
+profiling = False
+
+if profiling:
+    import resource
+
+##################################################################
+
+version = "2.99"
+rsmc_banner = """
+Reaction Systems SMT-Based Model Checking
+
+Version:  """ + version + """
+Author:   Artur Meski <artur.meski@gmail.com>
+"""
+
+##################################################################
+
+
+def print_banner():
+    print()
+    for line in rsmc_banner.split("\n"):
+        print(colour_str(C_GREEN, " " + 3 * "-" + " "), line)
+    print()
+
+##################################################################
+
+
+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")
+    parser.add_argument(
+        "-n", "--scaling-parameter",
+        help="scaling parameter value (used in some benchmarks)")
+    parser.add_argument("-s", "--special_mode",
+                        help="special mode (used in some benchmarks)")
+
+    args = parser.parse_args()
+
+    print_banner()
+    rs_testing.run_tests(args)
+
+##################################################################
+
+
+if __name__ == "__main__":
+    try:
+        if profiling:
+            import profile
+            profile.run('main()')
+        else:
+            main()
+    except KeyboardInterrupt:
+        print("\nQuitting...")
+        sys.exit(99)
+
+# EOF

reactics-smt/scripts/avg.py

@@ -0,0 +1,69 @@
+#!/usr/bin/env python
+
+import sys
+from os import listdir
+
+
+def proc_files(files):
+
+    fs = []
+    for f in files:
+        try:
+            fs.append(open(f, "r"))
+        except Exception as exct:
+            print("Failed to read file {:s}".format(f))
+            sys.exit(1)
+
+    done = False
+    result = []
+    while not done:
+
+        index = "NONE"
+        val_sum = 0
+
+        for f in fs:
+            try:
+                line = f.__next__()
+            except Exception:
+                done = True
+                break
+
+            sline = line.split()
+
+            assert(len(sline) == 2)
+
+            index, value = sline
+            value = float(value)
+
+            val_sum += value
+
+            # print(sline)
+
+        else:
+            val_avg = val_sum/len(fs)
+            result.append("{:s}\t{:f}".format(index, val_avg))
+
+    result.append("")
+    return result
+
+
+def proc_dirs(dirs):
+
+    first_dir = dirs[0]
+
+    for f in listdir(first_dir):
+        avg_out = proc_files(["{:s}/{:s}".format(d, f) for d in dirs])
+
+        with open(f, "w") as outfile:
+            outfile.write("\n".join(avg_out))
+
+
+def main():
+
+    proc_dirs(sys.argv[1:])
+
+    # proc_files(sys.argv[1:])
+
+
+if __name__ == "__main__":
+    main()

reactics-smt/scripts/mkplotdat.sh

@@ -0,0 +1,8 @@
+#!/bin/sh
+if [ "$1" = "" ];then
+	echo "Usage: $0 <input file>"
+	exit 1
+fi
+f="$1"
+cat $f | tr -d '()' | tr ',' '\t'
+

reactics-smt/scripts/run.sh

@@ -0,0 +1,7 @@
+#!/bin/sh
+
+for i in 2 3 4 5 6 7 8 9 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100;
+do
+  ./rssmt.py $i | tee "out/result_$i.out" 
+  echo "$i: $(tail -1 out/result_$i.out)" >> out/all.out
+done

reactics-smt/scripts/run_3d.sh

@@ -0,0 +1,24 @@
+#!/bin/sh
+
+# sudo systemsetup -getcomputersleep 
+sudo systemsetup -setcomputersleep Never
+# sudo systemsetup -setcomputersleep 1
+
+
+
+for i in `seq 2 2 20`
+do
+  #./rssmt.py $i | tee "out/result_$i.out" 
+  #echo "$i: $(tail -1 out/result_$i.out)" >> out/all.out
+  echo $i
+  for j in `seq 2 2 20`
+  do
+    echo $i $j
+	for t in 1 0
+	do
+	    ./rssmt.py $i $j $t
+	done
+  done
+done
+
+

reactics-smt/scripts/run_mutex_param.sh

@@ -0,0 +1,21 @@
+#!/bin/sh
+
+# sudo systemsetup -getcomputersleep
+#sudo systemsetup -setcomputersleep Never
+# sudo systemsetup -setcomputersleep 1
+
+for i in `seq 2 50`
+do
+    for special_mode in 1 2 3
+    do
+        echo "$i (sm=${special_mode})"
+        if [[ $special_mode -eq 1 ]]
+        then
+            ./rssmt.py -n $i -s $special_mode
+            ./rssmt.py -o -n $i -s $special_mode
+        else
+            ./rssmt.py -n $i -s $special_mode
+        fi
+    done
+done
+

reactics-smt/scripts/run_scalable_chain_rsLTL_FI.sh

@@ -0,0 +1,34 @@
+#!/bin/sh
+
+# sudo systemsetup -getcomputersleep
+sudo systemsetup -setcomputersleep Never
+# sudo systemsetup -setcomputersleep 1
+
+formula=1
+for i in `seq 2 50`
+do
+	echo $i
+	./rssmt.py $i 2 $formula
+done
+
+formula=2
+for i in `seq 2 50`
+do
+	echo $i
+	./rssmt.py $i 2 $formula
+done
+
+formula=3
+for i in `seq 2 50`
+do
+	echo $i
+	./rssmt.py $i 2 $formula
+done
+
+formula=4
+for i in `seq 2 50`
+do
+	echo $i
+	./rssmt.py $i 2 $formula
+done
+

reactics-smt/smt/__init__.py

@@ -0,0 +1,4 @@
+#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

reactics-smt/smt/smt_checker_rs.py

@@ -0,0 +1,314 @@
+"""
+SMT-based Model Checking Module for RS with Context Automaton
+"""
+
+from z3 import *
+from time import time
+from sys import stdout
+import resource
+
+
+class SmtCheckerRS(object):
+
+    def __init__(self, rsca):
+
+        rsca.sanity_check()
+
+        self.rs = rsca.rs
+        self.ca = rsca.ca
+
+        self.v = []
+        self.v_ctx = []
+        self.ca_state = []
+        self.next_level_to_encode = 0
+
+        self.solver = Solver()
+
+        self.verification_time = None
+
+    def prepare_all_variables(self):
+        """Encodes all the variables"""
+
+        self.prepare_state_variables()
+        self.prepare_context_variables()
+        self.next_level_to_encode += 1
+
+    def prepare_context_variables(self):
+        """Encodes all the context variables"""
+
+        level = self.next_level_to_encode
+
+        variables = []
+        for entity in self.rs.background_set:
+            variables.append(Bool("C"+str(level)+"_"+entity))
+
+        self.v_ctx.append(variables)
+
+    def prepare_rs_state_variables(self):
+        """Encodes all the state variables of the reaction system"""
+
+        level = self.next_level_to_encode
+
+        variables = []
+        for entity in self.rs.background_set:
+            variables.append(Bool("L"+str(level)+"_"+entity))
+        self.v.append(variables)
+
+    def prepare_context_controller_variables(self):
+        """Encodes all the variables required for controlling context sequences"""
+
+        level = self.next_level_to_encode
+
+        self.ca_state.append(Int("CA"+str(level)+"_state"))
+
+    def prepare_state_variables(self):
+        """Encodes all the state variables"""
+
+        self.prepare_rs_state_variables()
+        self.prepare_context_controller_variables()
+
+    def enc_rs_init_state(self, level):
+        """Encodes the initial state for the reaction system"""
+
+        rs_init_state_enc = True
+        for v in self.v[level]:
+            # the initial state is empty
+            rs_init_state_enc = simplify(And(rs_init_state_enc, Not(v)))
+        return rs_init_state_enc
+
+    def enc_context_controller_init_state(self, level):
+        """Encodes the initial state for controlling context sequences"""
+
+        return self.ca_state[level] == self.ca.get_init_state_id()
+
+    def enc_init_state(self, level):
+        """Encodes the initial state at the given level"""
+
+        init_state_enc = simplify(And(self.enc_rs_init_state(
+            level), self.enc_context_controller_init_state(level)))
+
+        return init_state_enc
+
+    def enc_enabledness(self, level, prod_entity):
+        """Encodes the enabledness condition for a given level and a given entity"""
+
+        rcts_for_prod_entity = self.rs.get_reactions_by_product()[prod_entity]
+
+        if rcts_for_prod_entity == []:
+            return False
+
+        enc_rct_prod = False
+        for reactants, inhibitors in rcts_for_prod_entity:
+            enc_reactants = True
+            enc_inhibitors = True
+            for reactant in reactants:
+                enc_reactants = simplify(And(enc_reactants, Or(
+                    self.v[level][reactant], self.v_ctx[level][reactant])))
+            for inhibitor in inhibitors:
+                enc_inhibitors = simplify(And(enc_inhibitors, Not(
+                    Or(self.v[level][inhibitor], self.v_ctx[level][inhibitor]))))
+
+            enc_rct_prod = simplify(
+                Or(enc_rct_prod, And(enc_reactants, enc_inhibitors)))
+
+        return enc_rct_prod
+
+    def enc_entity_production(self, level, prod_entity):
+        """Encodes the production of a given entity from a given level at level+1"""
+
+        enc_enab_cond = self.enc_enabledness(level, prod_entity)
+
+        enc_ent_prod = Or(
+            And(enc_enab_cond, self.v[level + 1][prod_entity]),
+            And(Not(enc_enab_cond),
+                Not(self.v[level + 1][prod_entity])))
+
+        return simplify(enc_ent_prod)
+
+    def enc_transition_relation(self, level):
+        """Encodes the combined transition relation"""
+
+        return simplify(
+            And(self.enc_rs_trans(level),
+                self.enc_automaton_trans(level)))
+
+    def enc_rs_trans(self, level):
+        """Encodes the transition relation"""
+
+        unused_entities = list(range(len(self.rs.background_set)))
+
+        enc_trans = True
+
+        for prod_entity in self.rs.get_reactions_by_product():
+            unused_entities.remove(prod_entity)
+
+            enc_trans = simplify(
+                And(enc_trans, self.enc_entity_production(level, prod_entity)))
+
+        for prod_entity in unused_entities:
+            enc_trans = simplify(
+                And(enc_trans, Not(self.v[level+1][prod_entity])))
+
+        return enc_trans
+
+    def enc_automaton_single_trans(self, level, transition):
+
+        src, ctx, dst = transition
+
+        src_enc = self.ca_state[level] == src
+        dst_enc = self.ca_state[level+1] == dst
+
+        all_ent = set(range(len(self.rs.background_set)))
+        incl_ctx = ctx
+        excl_ctx = all_ent - incl_ctx
+
+        ctx_enc = True
+
+        for c in incl_ctx:
+            ctx_enc = simplify(And(ctx_enc, self.v_ctx[level][c]))
+        for c in excl_ctx:
+            ctx_enc = simplify(And(ctx_enc, Not(self.v_ctx[level][c])))
+
+        enc_single_trans = simplify(And(src_enc, ctx_enc, dst_enc))
+
+        return enc_single_trans
+
+    def enc_automaton_trans(self, level):
+        """Encodes the transition relation for the context automaton"""
+
+        enc_trans = False
+        for transition in self.ca.transitions:
+            enc_trans = simplify(
+                Or(enc_trans, self.enc_automaton_single_trans(level, transition)))
+
+        return enc_trans
+
+    def enc_state(self, level, state):
+        """Encodes the state at the given level"""
+
+        enc = True
+
+        state_ids = self.rs.get_state_ids(state)
+
+        for entity in state_ids:
+            enc = And(enc, self.v[level][entity])
+
+        not_in_state = set(range(len(self.rs.background_set)))
+        not_in_state = not_in_state.difference(set(state_ids))
+
+        for entity in not_in_state:
+            enc = And(enc, Not(self.v[level][entity]))
+
+        return enc
+
+    def enc_non_exclusive_state(self, level, state):
+        """Encodes the state at the given level"""
+
+        enc = True
+
+        state_ids = self.rs.get_state_ids(state)
+
+        for entity in state_ids:
+            enc = And(enc, self.v[level][entity])
+
+        return 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
+
+        required_ids = self.rs.get_state_ids(required)
+        blocked_ids = self.rs.get_state_ids(blocked)
+
+        for e in required_ids:
+            enc = And(enc, self.v[level][e])
+        for e in blocked_ids:
+            enc = And(enc, Not(self.v[level][e]))
+
+        return simplify(enc)
+
+    def decode_witness(self, max_level, print_model=False):
+
+        m = self.solver.model()
+
+        if print_model:
+            print(m)
+
+        for level in range(max_level+1):
+
+            print("\n[Level=" + repr(level) + "]")
+
+            print("  State: {", end=""),
+            for var_id in range(len(self.v[level])):
+                if repr(m[self.v[level][var_id]]) == "True":
+                    print(" " + self.rs.get_entity_name(var_id), end="")
+            print(" }")
+
+            if level != max_level:
+                print("  Context set: ", end="")
+                print("{", end="")
+                for var_id in range(len(self.v[level])):
+                    if repr(m[self.v_ctx[level][var_id]]) == "True":
+                        print(" " + self.rs.get_entity_name(var_id), end="")
+                print(" }")
+
+    def check_reachability(
+            self, state, print_witness=True, print_time=True, print_mem=True):
+        """Main testing function"""
+
+        if not type(state) is tuple:
+            state = (state, [])
+
+        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))
+        current_level = 0
+
+        while True:
+            print("-----[ Working at level=" + str(current_level) + " ]-----")
+            stdout.flush()
+
+            self.prepare_all_variables()
+
+            # reachability test:
+            print("[i] Adding the reachability test...")
+            self.solver.push()
+            self.solver.add(self.enc_state_with_blocking(current_level, state))
+
+            result = self.solver.check()
+            if result == sat:
+                print("\n[+] SAT at level=" + str(current_level))
+                if print_witness:
+                    self.decode_witness(current_level)
+                break
+            else:
+                self.solver.pop()
+
+            print("[i] Unrolling the transition relation")
+            self.solver.add(self.enc_transition_relation(current_level))
+
+            print("-----[ level=" + str(current_level) + " done ]")
+            current_level += 1
+
+        if print_time:
+            # stop = time()
+            stop = resource.getrusage(resource.RUSAGE_SELF).ru_utime
+            self.verification_time = stop-start
+            print()
+            print("[i] Time: " + repr(self.verification_time))
+
+        if print_mem:
+            print(
+                "[i] Memory: " +
+                repr(
+                    resource.getrusage(resource.RUSAGE_SELF).ru_maxrss /
+                    (1024 * 1024)) + " MB")
+
+    def get_verification_time(self):
+        return self.verification_time

reactics-smt/smt/smt_checker_rsc.py

@@ -0,0 +1,670 @@
+"""
+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
+
+# def simplify(x):
+#     return x
+
+
+class SmtCheckerRSC(object):
+
+    def __init__(self, rsca):
+
+        rsca.sanity_check()
+
+        if not rsca.is_with_concentrations():
+            raise RuntimeError("RS and CA with concentrations expected")
+
+        self.rs = rsca.rs
+        self.ca = rsca.ca
+
+        self.initialise()
+
+    def initialise(self):
+        """Initialises all the variables used by the checker"""
+
+        self.v = []
+        self.v_ctx = []
+        self.ca_state = []
+        self.next_level_to_encode = 0
+
+        self.loop_position = Int("loop_position")
+
+        self.solver = Solver()  # For("QF_FD")
+
+        self.verification_time = None
+
+    def reset(self):
+        """Reinitialises the state of the checker"""
+
+        self.initialise()
+
+    def prepare_all_variables(self):
+        """Encodes all the variables"""
+
+        self.prepare_state_variables()
+        self.prepare_context_variables()
+        self.next_level_to_encode += 1
+
+    def prepare_context_variables(self):
+        """Encodes all the context variables"""
+
+        level = self.next_level_to_encode
+
+        variables = []
+        for entity in self.rs.background_set:
+            variables.append(Int("C"+str(level)+"_"+entity))
+
+        self.v_ctx.append(variables)
+
+    def prepare_state_variables(self):
+        """Encodes all the state variables"""
+
+        level = self.next_level_to_encode
+
+        variables = []
+        for entity in self.rs.background_set:
+            variables.append(Int("L"+str(level)+"_"+entity))
+        self.v.append(variables)
+
+        self.ca_state.append(Int("CA"+str(level)+"_state"))
+
+    def enc_concentration_levels_assertion(self, level):
+        """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.
+        """
+
+        enc_nz = True
+
+        for e_i in range(len(self.rs.background_set)):
+            v = self.v[level][e_i]
+            v_ctx = self.v_ctx[level][e_i]
+            e_max = self.rs.get_max_concentration_level(e_i)
+            enc_nz = simplify(
+                And(enc_nz, v >= 0, v_ctx >= 0, v <= e_max, v_ctx <= e_max))
+
+        return enc_nz
+
+    def enc_init_state(self, level):
+        """Encodes the initial state at the given level"""
+
+        rs_init_state_enc = True
+
+        for v in self.v[level]:
+            # the initial concentration levels are zeroed
+            rs_init_state_enc = simplify(And(rs_init_state_enc, v == 0))
+
+        ca_init_state_enc = self.ca_state[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_produced_concentration(self, level, prod_entity):
+        """Encodes the produced concentrations for the given level and entity"""
+
+        rcts_for_prod_entity = []
+        if prod_entity in self.rs.get_reactions_by_product():
+            rcts_for_prod_entity = self.rs.get_reactions_by_product()[
+                prod_entity]
+
+        meta_reactions = []
+        if prod_entity in self.rs.meta_reactions:
+            meta_reactions = self.rs.meta_reactions[prod_entity]
+
+        permanency_inhibition = None
+        if prod_entity in self.rs.permanent_entities:
+            permanency_inhibition = self.rs.permanent_entities[prod_entity]
+
+        if rcts_for_prod_entity == [] and meta_reactions == []:
+            # this should never happen
+            return simplify(self.v[level+1][prod_entity] == 0)
+
+        enc_enabledness = False
+
+        # ----------- ordinary reactions --------------------------------------------
+
+        enc_rct_prod = False
+
+        enc_ordinary_reactions_enabledness = False
+
+        for reactants, inhibitors, products in rcts_for_prod_entity:
+
+            enc_reactants = True
+            for reactant, concentration in reactants:
+                enc_reactants = simplify(And(enc_reactants, Or(
+                    self.v[level][reactant] >= concentration, self.v_ctx[level][reactant] >= concentration)))
+
+            enc_inhibitors = True
+            for inhibitor, concentration in inhibitors:
+                enc_inhibitors = simplify(And(enc_inhibitors, And(
+                    self.v[level][inhibitor] < concentration, self.v_ctx[level][inhibitor] < concentration)))
+
+            enc_rct_enabled = And(enc_reactants, enc_inhibitors)
+            enc_products = self.v[level+1][products[0][0]] == products[0][1]
+            enc_rct_prod = simplify(
+                If(enc_rct_enabled, enc_products, enc_rct_prod))
+            enc_enabledness = simplify(Or(enc_enabledness, enc_rct_enabled))
+
+            enc_ordinary_reactions_enabledness = simplify(
+                Or(enc_ordinary_reactions_enabledness, enc_rct_enabled))
+
+        # -------- meta reactions ---------------------------------------------------
+
+        for r_type, command_entity, reactants, inhibitors in meta_reactions:
+
+            # command entity is e.g. 'inc' for incrementation operation
+            # (inc,W) gives us the value W by which the given entity's value should be incremented
+
+            enc_reactants = True
+            enc_inhibitors = True
+
+            for reactant, concentration in reactants:
+                enc_reactants = simplify(And(enc_reactants, Or(
+                    self.v[level][reactant] >= concentration, self.v_ctx[level][reactant] >= concentration)))
+
+            # command entity needs to be present (with concentration level > 0) in order to perform the operation
+            enc_reactants = simplify(And(enc_reactants, Or(
+                self.v[level][command_entity] > 0, self.v_ctx[level][command_entity] > 0)))
+
+            for inhibitor, concentration in inhibitors:
+                enc_inhibitors = simplify(And(enc_inhibitors, And(
+                    self.v[level][inhibitor] < concentration, self.v_ctx[level][inhibitor] < concentration)))
+
+            if r_type == "inc":
+                value_after_inc = If(
+                    self.v[level][prod_entity] > self.v_ctx[level]
+                    [prod_entity],
+                    self.v[level][prod_entity],
+                    self.v_ctx[level][prod_entity]) + If(
+                    self.v[level][command_entity] > self.
+                    v_ctx[level][command_entity],
+                    self.v[level][command_entity],
+                    self.v_ctx[level][command_entity])
+                enc_products = self.v[level+1][prod_entity] == value_after_inc
+
+            elif r_type == "dec":
+                value_after_dec = simplify(
+                    If(
+                        self.v[level][prod_entity] >
+                        self.v_ctx[level]
+                        [prod_entity],
+                        self.v[level][prod_entity],
+                        self.v_ctx[level]
+                        [prod_entity]) -
+                    If(
+                        self.v[level]
+                        [command_entity] > self.
+                        v_ctx[level][command_entity],
+                        self.v[level]
+                        [command_entity],
+                        self.v_ctx[level]
+                        [command_entity]))
+                enc_products = self.v[level+1][prod_entity] == If(
+                    value_after_dec < 0, 0, value_after_dec)
+
+            else:
+                raise RuntimeError(
+                    "Unknown meta-reaction type: " + repr(r_type))
+
+            enc_meta_reaction_enabledness = And(
+                enc_reactants, enc_inhibitors,
+                Not(enc_ordinary_reactions_enabledness))
+            enc_enabledness = simplify(
+                Or(enc_enabledness, enc_meta_reaction_enabledness))
+            enc_rct_prod = simplify(Or(enc_rct_prod, And(
+                enc_meta_reaction_enabledness, enc_products)))
+
+        # -----------------------------------------------------------------------------
+
+        if not permanency_inhibition == None:
+
+            enc_reactants = Or(self.v[level][prod_entity] >= concentration,
+                               self.v_ctx[level][prod_entity] >= concentration)
+
+            enc_inhibitors = True
+            for inhibitor, concentration in permanency_inhibition:
+                enc_inhibitors = simplify(And(enc_inhibitors, And(
+                    self.v[level][inhibitor] < concentration, self.v_ctx[level][inhibitor] < concentration)))
+            enc_products = simplify(
+                self.v[level + 1][prod_entity] ==
+                If(
+                    self.v[level][prod_entity] > self.v_ctx[level]
+                    [prod_entity],
+                    self.v[level][prod_entity],
+                    self.v_ctx[level][prod_entity]))
+
+            enc_permanency_enabledness = And(
+                enc_reactants, enc_inhibitors,
+                Not(enc_ordinary_reactions_enabledness))
+            enc_enabledness = simplify(
+                Or(enc_enabledness, enc_permanency_enabledness))
+            enc_permanency = And(enc_permanency_enabledness, enc_products)
+            enc_rct_prod = simplify(Or(enc_rct_prod, enc_permanency))
+
+        # -----------------------------------------------------------------------------
+
+        enc_when_to_produce_zero_conc = simplify(
+            And(Not(enc_enabledness), self.v[level+1][prod_entity] == 0))
+
+        enc_rct_prod = Or(enc_rct_prod, enc_when_to_produce_zero_conc)
+        return enc_rct_prod
+
+    def enc_transition_relation(self, level):
+        return simplify(
+            And(self.enc_rs_trans(level),
+                self.enc_automaton_trans(level)))
+
+    def enc_rs_trans(self, level):
+        """Encodes the transition relation"""
+
+        unused_entities = set(range(len(self.rs.background_set)))
+
+        enc_trans = True
+
+        reactions = self.rs.get_reactions_by_product()
+        meta_reactions = self.rs.meta_reactions
+
+        for prod_entity in chain(reactions, meta_reactions):
+            unused_entities.discard(prod_entity)
+            enc_trans = simplify(
+                And(enc_trans, self.enc_produced_concentration(level, prod_entity)))
+
+        for prod_entity in unused_entities:
+            enc_trans = simplify(
+                And(enc_trans, self.v[level+1][prod_entity] == 0))
+
+        return enc_trans
+
+    def enc_automaton_trans(self, level):
+        """Encodes the transition relation for the context automaton"""
+
+        enc_trans = False
+
+        for src, ctx, dst in self.ca.transitions:
+            src_enc = self.ca_state[level] == src
+            dst_enc = self.ca_state[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.v_ctx[level][e] == c))
+
+            for e in excl_ctx:
+                ctx_enc = simplify(And(ctx_enc, self.v_ctx[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)
+
+        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, print_model=False):
+
+        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.v[level])):
+                var_rep = repr(m[self.v[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.v[level])):
+                    var_rep = repr(m[self.v_ctx[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(" }")
+
+    def check_rsltl(
+            self, formula, print_witness=True, print_time=True, print_mem=True,
+            max_level=None):
+        """Bounded Model Checking for rsLTL properties"""
+
+        self.reset()
+
+        print("[" + colour_str(C_BOLD, "i") +
+              "] Running rsLTL bounded model checking")
+        print("[" + colour_str(C_BOLD, "i") + "] Formula: " + str(formula))
+
+        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))
+
+        encoder = rsLTL_Encoder(self)
+        encoder.load_variables(
+            var_rs=self.v,
+            var_ctx=self.v_ctx,
+            var_loop_pos=self.loop_position)
+
+        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:
+            self.solver.push()
+
+            print("[" + colour_str(C_BOLD, "i") +
+                  "] Generating the formula encoding...")
+
+            f = encoder.get_encoding(formula, self.current_level)
+            ncalls = encoder.get_ncalls()
+
+            print("[" + colour_str(C_BOLD, "i") + "] Cache hits: " +
+                  str(encoder.get_cache_hits()) + ", encode calls: " +
+                  str(ncalls[0]) + " (approx: " + str(ncalls[1]) + ")")
+            print("[" + colour_str(C_BOLD, "i") +
+                  "] Adding the formula to the solver...")
+
+            encoder.flush_cache()
+            self.solver.add(f)
+
+            print("[" + colour_str(C_BOLD, "i") +
+                  "] Adding the loops encoding...")
+            self.solver.add(self.get_loop_encodings())
+
+            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 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 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