diff --git a/colour.py b/colour.py index 07cc048..80b7c65 100644 --- a/colour.py +++ b/colour.py @@ -13,5 +13,16 @@ C_MARK_ERROR = C_BOLD + "[" + C_RED + "!" + C_ENDC + C_BOLD + "]" + C_ENDC def colour_str(col, s): return col + s + C_ENDC +def green_str(s): + return C_GREEN + s + C_ENDC + def print_error(s): - print(C_MARK_ERROR + " " + C_RED + s + C_ENDC) \ No newline at end of file + 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 diff --git a/logics/bags.py b/logics/bags.py index 30ef358..fd386e7 100644 --- a/logics/bags.py +++ b/logics/bags.py @@ -1,26 +1,29 @@ from enum import Enum -BagDesc_oper = Enum('BagDesc_oper', 'entity true l_and l_or l_not lt le eq ge gt') +BagDesc_oper = Enum( + 'BagDesc_oper', 'entity true l_and l_or l_not lt le eq ge gt') + class BagDescription(object): - def __init__(self, f_type, L_oper = None, R_oper = None, entity = ""): + + 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) + " )" + 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) + " )" + return "(" + repr(self.left_operand) + " | " + repr(self.right_operand) + ")" if self.f_type == BagDesc_oper.l_not: - return "~" + repr(self.left_operand) + return "~(" + repr(self.left_operand) + ")" if self.f_type == BagDesc_oper.lt: return repr(self.left_operand) + " < " + repr(self.right_operand) if self.f_type == BagDesc_oper.le: @@ -31,45 +34,56 @@ class BagDescription(object): 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: " + str(operand)) - + if not( + isinstance(operand, BagDescription) + or isinstance(operand, int)): + raise RuntimeError( + "Unexpected operand type for a bag: {:s} (type: {:s})".format( + str(operand), str(type(operand)))) + @classmethod def f_entity(cls, entity_name): - return cls(BagDesc_oper.entity, entity = 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) - + 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) - + 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) - + 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) - + 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) - + 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) - + 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) - + return BagDescription(BagDesc_oper.l_or, L_oper=self, R_oper=other) + def __invert__(self): - return BagDescription(BagDesc_oper.l_not, L_oper = self) - + return BagDescription(BagDesc_oper.l_not, L_oper=self) + # EOF diff --git a/logics/rsltl.py b/logics/rsltl.py index e8bd959..91fdcc1 100644 --- a/logics/rsltl.py +++ b/logics/rsltl.py @@ -3,43 +3,61 @@ from enum import Enum from logics.bags import * -rsLTL_form_type = Enum('rsLTL_form_type', 'bag l_and l_or l_not l_implies t_globally t_finally t_next t_until t_release') +rsLTL_form_type = Enum( + 'rsLTL_form_type', + 'bag l_and l_or l_not l_implies t_globally t_finally t_next t_until t_release') + class Formula_rsLTL(object): - - def __init__(self, f_type, L_oper = None, R_oper = None, sub_oper = None, bag = None): + + def __init__( + self, f_type, L_oper=None, R_oper=None, sub_oper=None, bag=None): self.f_type = f_type self.left_operand = L_oper self.right_operand = R_oper self.sub_operand = sub_oper self.bag_descr = bag - + + # it's silly, but it's a frequent mistake + if callable(sub_oper): + raise RuntimeError( + "sub_oper should not be a function. Missing () for {0}?".format(sub_oper)) + if isinstance(self.left_operand, BagDescription): self.left_operand = Formula_rsLTL.f_bag(self.left_operand) if isinstance(self.right_operand, BagDescription): self.right_operand = Formula_rsLTL.f_bag(self.right_operand) - + + if self.sub_operand is True: + self.sub_operand = BagDescription.f_TRUE() + def __repr__(self): if self.f_type == rsLTL_form_type.bag: return repr(self.bag_descr) if self.f_type == rsLTL_form_type.l_not: - return "~( " + repr(self.left_operand) + " )" + return "~(" + repr(self.left_operand) + ")" if self.f_type == rsLTL_form_type.t_globally: - return "G[" + repr(self.sub_operand) + "]( " + repr(self.left_operand) + " )" + 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) + " )" + 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) + " )" + 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) + " )" + 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) + " )" + 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) + " )" + return "(" + repr(self.left_operand) + " => " + repr(self.right_operand) + ")" if self.f_type == rsLTL_form_type.t_until: - return "( " + repr(self.left_operand) + " U[" + repr(self.sub_operand) + "] " + repr(self.right_operand) + " )" + return "(" + repr( + self.left_operand) + " U[" + repr( + self.sub_operand) + "] " + repr( + self.right_operand) + ")" if self.f_type == rsLTL_form_type.t_release: - return "( " + repr(self.left_operand) + " R[" + repr(self.sub_operand) + "] " + repr(self.right_operand) + " )" + return "(" + repr( + self.left_operand) + " R[" + repr( + self.sub_operand) + "] " + repr( + self.right_operand) + ")" @property def is_bag(self): @@ -47,48 +65,55 @@ class Formula_rsLTL(object): @classmethod def f_bag(cls, bag_descr): - return cls(rsLTL_form_type.bag, bag = 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) - + 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) - + 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) - + 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) + 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) - + 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) + 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) + 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) - + 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) - + 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) - + 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) + return Formula_rsLTL(rsLTL_form_type.l_not, L_oper=self) -# EOF \ No newline at end of file +# EOF diff --git a/logics/rsltl_encoder.py b/logics/rsltl_encoder.py index 79844b8..222535c 100644 --- a/logics/rsltl_encoder.py +++ b/logics/rsltl_encoder.py @@ -8,16 +8,34 @@ class rsLTL_Encoder(object): def __init__(self, smt_checker): self.smt_checker = smt_checker - self.v = smt_checker.v - self.v_ctx = smt_checker.v_ctx self.rs = smt_checker.rs - self.loop_position = smt_checker.loop_position + + self.v = None + self.v_ctx = None + self.loop_position = None + + # self.load_variables( + # var_rs=smt_checker.v, + # var_ctx=smt_checker.v_ctx, + # var_loop_pos=smt_checker.loop_position) self.init_ncalls() + + def load_variables(self, var_rs, var_ctx, var_loop_pos): + + self.v = var_rs + self.v_ctx = var_ctx + self.loop_position = var_loop_pos def get_encoding(self, formula, bound): + + assert self.v is not None + assert self.v_ctx is not None + assert self.loop_position is not None + self.cache_init(bound) self.init_ncalls() + return self.encode(formula, 0, bound) def init_ncalls(self): @@ -102,12 +120,18 @@ class rsLTL_Encoder(object): if bag_formula.f_type == BagDesc_oper.gt: return self.encode_bag(bag_formula.left_operand, level, context) > int(bag_formula.right_operand) + + assert False, "Unsupported case {:s}".format(bag_formula.f_type) def encode_bag_state(self, bag_formula, level): - return self.encode_bag(bag_formula, level) + res = self.encode_bag(bag_formula, level) + assert res is not None + return res def encode_bag_ctx(self, bag_formula, level): - return self.encode_bag(bag_formula, level, context=True) + res = self.encode_bag(bag_formula, level, context=True) + assert res is not None + return res def encode(self, formula, level, bound): @@ -149,7 +173,7 @@ class rsLTL_Encoder(object): elif formula.f_type == rsLTL_form_type.l_implies: enc = Implies( - self.encode(formula.left_operand, level, bound), + self.encode(formula.left_operand, level, bound), self.encode(formula.right_operand, level, bound) ) diff --git a/rs/__init__.py b/rs/__init__.py index b500f08..129a278 100644 --- a/rs/__init__.py +++ b/rs/__init__.py @@ -2,10 +2,11 @@ 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 \ No newline at end of file +from rs.reaction_system_with_autnet import ReactionSystemWithNetworkOfAutomata diff --git a/rs/reaction_system_with_automaton.py b/rs/reaction_system_with_automaton.py index 9690ce2..500a883 100644 --- a/rs/reaction_system_with_automaton.py +++ b/rs/reaction_system_with_automaton.py @@ -1,4 +1,5 @@ from rs.reaction_system_with_concentrations import ReactionSystemWithConcentrations +from rs.reaction_system_with_concentrations_param import ReactionSystemWithConcentrationsParam from rs.context_automaton_with_concentrations import ContextAutomatonWithConcentrations class ReactionSystemWithAutomaton(object): @@ -11,6 +12,17 @@ class ReactionSystemWithAutomaton(object): self.rs.show(soft) self.ca.show() + def is_concentr_and_param_compatible(self): + """ + Checks if the underlying RS/CA are compatible + with parameters and concentrations + """ + if not isinstance(self.rs, ReactionSystemWithConcentrationsParam): + return False + if not isinstance(self.ca, ContextAutomatonWithConcentrations): + return False + return True + def is_with_concentrations(self): if not isinstance(self.rs, ReactionSystemWithConcentrations): return False diff --git a/rs/reaction_system_with_concentrations.py b/rs/reaction_system_with_concentrations.py index 4167aa2..daf9fde 100644 --- a/rs/reaction_system_with_concentrations.py +++ b/rs/reaction_system_with_concentrations.py @@ -56,8 +56,7 @@ class ReactionSystemWithConcentrations(ReactionSystem): if len(e) == 2 and type(e[1]) is int: return True - print("FATAL. Invalid entity+concentration:") - print(e) + print("FATAL. Invalid entity+concentration: {:s}".format(e)) exit(1) return False diff --git a/rs/reaction_system_with_concentrations_param.py b/rs/reaction_system_with_concentrations_param.py new file mode 100644 index 0000000..4ed6726 --- /dev/null +++ b/rs/reaction_system_with_concentrations_param.py @@ -0,0 +1,529 @@ +from sys import exit +from colour import * + +from rs.reaction_system import ReactionSystem + +class ParameterObj(object): + + def __init__(self, name): + self.name = name + + def __repr__(self): + return "@{0}".format(self.name) + +def is_param(some_object): + if isinstance(some_object, ParameterObj): + return True + else: + return False + +class ReactionSystemWithConcentrationsParam(ReactionSystem): + + def __init__(self): + + self.reactions = [] + self.parameters = dict() + self.meta_reactions = dict() + self.permanent_entities = dict() + self.background_set = [] + self.context_entities = [] # legacy. to be removed + self.reactions_by_prod = None + self.max_concentration = 0 + self.max_conc_per_ent = dict() + + def add_bg_set_entity(self, e): + name = "" + def_max_conc = -1 + if type(e) is tuple and len(e) == 2: + name, def_max_conc = e + elif type(e) is str: + name = e + print("\nWARNING: no maximal concentration level specified for:", e, "\n") + else: + raise RuntimeError( + "Bad entity type when adding background set element") + + self.assume_not_in_bgset(name) + self.background_set.append(name) + + if def_max_conc != -1: + ent_id = self.get_entity_id(name) + self.max_conc_per_ent.setdefault(ent_id, 0) + if self.max_conc_per_ent[ent_id] < def_max_conc: + self.max_conc_per_ent[ent_id] = def_max_conc + if self.max_concentration < def_max_conc: + self.max_concentration = def_max_conc + + def get_param(self, name): + if self.has_param(name): + return self.parameters[name] + else: + param = ParameterObj(name) + self.add_param(param) + return param + + def has_param(self, name): + return name in self.parameters + + def add_param(self, param): + if param in self.parameters: + raise RuntimeError("Parameter {:s} already exists".format(param)) + param_key = param.name + self.parameters[param_key] = param + self.parameters[param_key].idx = len(self.parameters) + + def get_max_concentration_level(self, e): + + if e in self.max_conc_per_ent: + return self.max_conc_per_ent[e] + else: + return self.max_concentration + + def is_valid_entity_with_concentration(self, e): + """Sanity check for entities with concentration""" + + if type(e) is tuple: + if len(e) == 2 and type(e[1]) is int: + return True + + if type(e) is list: + if len(e) == 2 and type(e[1]) is int: + return True + + print("FATAL. Invalid entity+concentration: {:s}".format(e)) + exit(1) + + return False + + def get_state_ids(self, state): + """Returns entities of the given state without levels""" + # return [e for e,c in state] + return [self.get_entity_id(e) for e in state] + + def has_non_zero_concentration(self, elem): + if elem[1] < 1: + raise RuntimeError( + "Unexpected concentration level in state: " + str(elem)) + + def process_rip(self, R, I, P, ignore_empty_R=False): + """Chcecks concentration levels and converts entities names into their ids""" + + if R == [] and not ignore_empty_R: + raise RuntimeError("No reactants defined") + + # + # REACTANTS + # + reactants = [] + if isinstance(R, ParameterObj): + reactants = R + else: + for r in R: + self.is_valid_entity_with_concentration(r) + self.has_non_zero_concentration(r) + entity, level = r + reactants.append((self.get_entity_id(entity), level)) + if self.max_concentration < level: + self.max_concentration = level + + # + # INHIBITORS + # + inhibitors = [] + if isinstance(I, ParameterObj): + inhibitors = I + else: + for i in I: + self.is_valid_entity_with_concentration(i) + self.has_non_zero_concentration(i) + entity, level = i + inhibitors.append((self.get_entity_id(entity), level)) + if self.max_concentration < level: + self.max_concentration = level + + # + # PRODUCTS + # + products = [] + if isinstance(P, ParameterObj): + products = P + else: + for p in P: + self.is_valid_entity_with_concentration(p) + self.has_non_zero_concentration(p) + entity, level = p + products.append((self.get_entity_id(entity), level)) + + return reactants, inhibitors, products + + def is_parametric_reaction(self, reaction): + result = any([isinstance(r_set, ParameterObj) for r_set in reaction]) + return result + + def add_reaction(self, R, I, P): + """Adds a reaction + + R, I, and P are sets of entities (not their IDs) + """ + + if P == []: + raise RuntimeError("No products defined") + reaction = self.process_rip(R,I,P) + + self.reactions.append(reaction) + + def add_reaction_without_reactants(self, R, I, P): + """Adds a reaction""" + + if P == []: + raise RuntimeError("No products defined") + reaction = self.process_rip(R,I,P,ignore_empty_R=True) + self.reactions.append(reaction) + + def add_reaction_inc(self, incr_entity, incrementer, R, I): + """Adds a macro/meta reaction for increasing the value of incr_entity""" + + reactants,inhibitors,products = self.process_rip(R,I,[],ignore_empty_R=True) + incr_entity_id = self.get_entity_id(incr_entity) + self.meta_reactions.setdefault(incr_entity_id,[]) + self.meta_reactions[incr_entity_id].append(("inc", self.get_entity_id(incrementer), reactants, inhibitors)) + + def add_reaction_dec(self, decr_entity, decrementer, R, I): + """Adds a macro/meta reaction for decreasing the value of incr_entity""" + + reactants,inhibitors,products = self.process_rip(R,I,[],ignore_empty_R=True) + decr_entity_id = self.get_entity_id(decr_entity) + self.meta_reactions.setdefault(decr_entity_id,[]) + self.meta_reactions[decr_entity_id].append(("dec", self.get_entity_id(decrementer), reactants, inhibitors)) + + def add_permanency(self, ent, I): + """Sets entity to be permanent unless it is inhibited""" + + ent_id = self.get_entity_id(ent) + + if ent_id in self.permanent_entities: + raise RuntimeError("Permanency for {0} already defined.".format(ent)) + + inhibitors = self.process_rip([],I,[],ignore_empty_R=True)[1] + self.permanent_entities[ent_id] = inhibitors + + def set_context_entities(self, entities): + raise NotImplementedError + + def entities_names_set_to_str(self, entities): + s = "" + for entity in entities: + s += entity + ", " + s = s[:-2] + return s + + def entities_ids_set_to_str(self, entities): + s = "" + for entity in entities: + s += self.get_entity_name(entity) + ", " + s = s[:-2] + return s + + def state_to_str(self, state): + """ + If state is a parameter, we return + the string representation of the whole state + which should be the name of the parameter + """ + if isinstance(state, ParameterObj): + return str(state) + else: + s = "" + for ent,level in state: + s += self.get_entity_name(ent) + "=" + str(level) + ", " + s = s[:-2] + return s + + def show_background_set(self): + print(C_MARK_INFO + " Background set: {" + self.entities_names_set_to_str(self.background_set) + "}") + + def show_meta_reactions(self): + print(C_MARK_INFO + " Meta reactions:") + for param_ent,reactions in self.meta_reactions.items(): + for r_type,command,reactants,inhibitors in reactions: + if r_type == "inc" or r_type == "dec": + print(" - [ Type=" + repr(r_type) + " Operand=( " + self.get_entity_name(param_ent) + \ + " ) Command=( " + self.get_entity_name(command) + " ) ] -- ( R={" + self.state_to_str(reactants) + "}, I={" + self.state_to_str(inhibitors) + "} )") + else: + raise RuntimeError("Unknown meta-reaction type: " + repr(r_type)) + + def show_max_concentrations(self): + print(C_MARK_INFO + " Maximal allowed concentration levels:") + for e,max_conc in self.max_conc_per_ent.items(): + print(" - {0:^20} = {1:<6}".format(self.get_entity_name(e),max_conc)) + + def show_permanent_entities(self): + print(C_MARK_INFO + " Permanent entities:") + for e,inhibitors in self.permanent_entities.items(): + print(" - {0:^20}{1:<6}".format(self.get_entity_name(e) + ": ","I={" + self.state_to_str(inhibitors) + "}")) + + def show(self, soft=False): + self.show_background_set() + self.show_reactions(soft) + # self.show_param_reactions(soft) + self.show_permanent_entities() + self.show_meta_reactions() + self.show_max_concentrations() + + def get_producible_entities(self): + """ + Returns the set of entities that appear as products of + reactions. + """ + + producible_entities = set() + + for reaction in self.reactions: + product_entities = [e for e,c in reaction[2] if c > 0] + producible_entities = producible_entities.union(set(product_entities)) + + return producible_entities + + def get_reactions_by_product(self): + """Sorts reactions by their products and returns a dictionary of products""" + + # assert False + + if self.reactions_by_prod != None: + return self.reactions_by_prod + + producible_entities = self.get_producible_entities() + + reactions_by_prod = {} + + for p_e in producible_entities: + reactions_by_prod[p_e] = [] + rcts_for_p_e = reactions_by_prod[p_e] + + for r in self.reactions: + product_entities = [e for e,c in r[2]] + + if p_e in product_entities: + reactants = r[0] + inhibitors = r[1] + products = [(e,c) for e,c in r[2] if e == p_e] + + prod_conc = products[0][1] + insert_place = None + + # we need to order the reactions w.r.t. the concentration levels produced (increasing order) + for i in range(0,len(rcts_for_p_e)): + + checked_conc = rcts_for_p_e[i][2][0][1] + if prod_conc <= checked_conc: + insert_place = i + break + + if insert_place == None: # empty or the is only one element which is smaller than the element being added + rcts_for_p_e.append((reactants, inhibitors, products)) # we append (to the end) + else: + rcts_for_p_e.insert(insert_place,(reactants, inhibitors, products)) + + + # save in cache + self.reactions_by_prod = reactions_by_prod + + return reactions_by_prod + + def get_reaction_system(self): + """ + Translates RSC into RS + """ + + rs = ReactionSystem() + + for reactants, inhibitors, products in self.reactions: + + new_reactants = [] + new_inhibitors = [] + new_products = [] + + for ent, conc in reactants: + n = self.get_entity_name(ent) + "#" + str(conc) + rs.ensure_bg_set_entity(n) + new_reactants.append(n) + + for ent, conc in inhibitors: + n = self.get_entity_name(ent) + "#" + str(conc) + rs.ensure_bg_set_entity(n) + new_inhibitors.append(n) + + for ent, conc in products: + for i in range(1, conc + 1): + n = self.get_entity_name(ent) + "#" + str(i) + rs.ensure_bg_set_entity(n) + new_products.append(n) + + rs.add_reaction(new_reactants, new_inhibitors, new_products) + + for param_ent, reactions in self.meta_reactions.items(): + for r_type, command, reactants, inhibitors in reactions: + + param_ent_name = self.get_entity_name(param_ent) + + new_reactants = [] + new_inhibitors = [] + + for ent, conc in reactants: + n = self.get_entity_name(ent) + "#" + str(conc) + rs.ensure_bg_set_entity(n) + new_reactants.append(n) + + for ent, conc in inhibitors: + n = self.get_entity_name(ent) + "#" + str(conc) + rs.ensure_bg_set_entity(n) + new_inhibitors.append(n) + + max_cmd_c = self.max_concentration + if command in self.max_conc_per_ent: + max_cmd_c = self.max_conc_per_ent[command] + else: + print( + "WARNING:\n\tThere is no maximal concentration level defined for " + + self.get_entity_name(command)) + print("\tThis is a very bad idea -- expect degraded performance\n") + + for l in range(1, max_cmd_c + 1): + + cmd_ent = self.get_entity_name(command) + "#" + str(l) + rs.ensure_bg_set_entity(cmd_ent) + + if r_type == "inc": + + # pre_conc -- predecessor concentration + # succ_conc -- successor concentration concentration + + for i in range(1, self.max_concentration): + pre_conc = param_ent_name + "#" + str(i) + rs.ensure_bg_set_entity(pre_conc) + new_products = [] + succ_value = i + l + for j in range(1, succ_value + 1): + if j > self.max_concentration: + break + new_p = param_ent_name + "#" + str(j) + rs.ensure_bg_set_entity(new_p) + new_products.append(new_p) + if new_products != []: + rs.add_reaction( + set(new_reactants + [pre_conc, cmd_ent]), + set(new_inhibitors), + set(new_products)) + + elif r_type == "dec": + for i in range(1, self.max_concentration + 1): + pre_conc = param_ent_name + "#" + str(i) + rs.ensure_bg_set_entity(pre_conc) + new_products = [] + succ_value = i - l + for j in range(1, succ_value + 1): + if j > self.max_concentration: + break + new_p = param_ent_name + "#" + str(j) + rs.ensure_bg_set_entity(new_p) + new_products.append(new_p) + if new_products != []: + rs.add_reaction( + set(new_reactants + [pre_conc, cmd_ent]), + set(new_inhibitors), + set(new_products)) + + else: + raise RuntimeError( + "Unknown meta-reaction type: " + repr(r_type)) + + for ent, inhibitors in self.permanent_entities.items(): + + max_c = self.max_concentration + if ent in self.max_conc_per_ent: + max_c = self.max_conc_per_ent[ent] + else: + print( + "WARNING:\n\tThere is no maximal concentration level defined for " + + self.get_entity_name(ent)) + print("\tThis is a very bad idea -- expect degraded performance\n") + + def e_value(i): + return self.get_entity_name(ent) + "#" + str(i) + + for value in range(1, max_c + 1): + + new_reactants = [] + new_inhibitors = [] + new_products = [] + + new_reactants = [e_value(value)] + + for e_inh, conc in inhibitors: + n = self.get_entity_name(e_inh) + "#" + str(conc) + rs.ensure_bg_set_entity(n) + new_inhibitors.append(n) + + for i in range(1, value + 1): + new_products.append(e_value(i)) + + rs.add_reaction(new_reactants, new_inhibitors, new_products) + + return rs + + +# class ReactionSystemWithAutomaton(object): +# +# def __init__(self, reaction_system, context_automaton): +# self.rs = reaction_system +# self.ca = context_automaton +# +# def show(self, soft=False): +# self.rs.show(soft) +# self.ca.show() +# +# def is_concentr_and_param_compatible(self): +# """ +# Checks if the underlying RS/CA are compatible +# with parameters and concentrations +# """ +# if not isinstance(self.rs, ReactionSystemWithConcentrationsParam): +# return False +# if not isinstance(self.ca, ContextAutomatonWithConcentrations): +# return False +# return True +# +# def is_with_concentrations(self): +# """ +# Checks if the underlying RS and CA provide +# concentration levels +# """ +# if not isinstance(self.rs, ReactionSystemWithConcentrations): +# return False +# if not isinstance(self.ca, ContextAutomatonWithConcentrations): +# return False +# return True +# +# def sanity_check(self): +# pass +# +# def get_ordinary_reaction_system_with_automaton(self): +# +# if not self.is_with_concentrations(): +# raise RuntimeError("Not RS/CA with concentrations") +# +# ors = self.rs.get_reaction_system() +# oca = self.ca.get_automaton_with_flat_contexts(ors) +# +# return ReactionSystemWithAutomaton(ors, oca) +# + +# class ReactionSystemWithConcentrationWithAutomaton(ReactionSystemWithAutomaton): +# +# def __init__(self, reaction_system, context_automaton): +# self.rs = reaction_system +# self.ca = context_automaton +# +# def show(self, soft=False): +# self.rs.show(soft) +# self.ca.show() diff --git a/rs_examples.py b/rs_examples.py index bbe9478..f0e26b9 100755 --- a/rs_examples.py +++ b/rs_examples.py @@ -187,3 +187,4 @@ def heat_shock_response(print_system=True,verify_rsc=True): def state_translate_rsc2rs(p): return [e[0] + "#" + str(e[1]) for e in p] +# EOF \ No newline at end of file diff --git a/rs_testing.py b/rs_testing.py index 6ecb83c..5d0e34b 100644 --- a/rs_testing.py +++ b/rs_testing.py @@ -2,81 +2,322 @@ from rs import * from smt import * import rs_examples from logics import * +from rsltl_shortcuts import * + import sys import resource -def run_tests(): +def run_tests(cmd_args): # test_extended_automaton() # process() # heat_shock_response() # scalable_chain(print_system=True) - example44() + # example44() + # example44_param() + # heat_shock_response_param(cmd_args) + # simple_param(cmd_args) + + gene_expression(cmd_args) + +def gene_expression(cmd_args): + """ + Simple gene expression example + """ + r = ReactionSystemWithConcentrationsParam() -def example44(): + r.add_bg_set_entity(("x", 1)) + r.add_bg_set_entity(("xp", 1)) + r.add_bg_set_entity(("X", 1)) + + r.add_bg_set_entity(("y", 1)) + r.add_bg_set_entity(("yp", 1)) + r.add_bg_set_entity(("Y", 1)) - r = ReactionSystemWithConcentrations() - r.add_bg_set_entity(("x",2)) - r.add_bg_set_entity(("y",4)) - r.add_bg_set_entity(("h",2)) - r.add_bg_set_entity(("m",1)) + r.add_bg_set_entity(("z", 1)) + r.add_bg_set_entity(("zp", 1)) + r.add_bg_set_entity(("Z", 1)) - r.add_reaction([("y",1),("x",1)], [("y",2),("h",1)], [("y",2)]) - r.add_reaction([("y",2),("x",2)], [("y",3),("h",1)], [("y",3)]) - r.add_reaction([("y",3),("h",1)], [("y",4),("h",2)], [("y",4)]) - r.add_reaction([("y",4),("h",1)], [("h",2)], [("y",3)]) - r.add_reaction([("y",4),("x",2)], [("h",1)], [("y",2)]) - r.add_reaction([("m",1)], [("y",3)], [("m",1)]) + r.add_bg_set_entity(("h", 1)) + r.add_bg_set_entity(("Q", 1)) + r.add_bg_set_entity(("U", 1)) + + all_entities = set(r.background_set) + + r.add_reaction([("x",1)],[("h",1)],[("x",1)]) + + # r.add_reaction([("x",1)],[("h",1)],[("xp",1)]) + param_p1 = r.get_param("P1") + r.add_reaction(param_p1,[("h",1)],[("xp",1)]) + + # r.add_reaction([("x",1),("xp",1)],[("h",1)],[("X",1)]) + param_p2_r = r.get_param("P2_r") + param_p2_i = r.get_param("P2_i") + param_p2_p = r.get_param("P2_p") + r.add_reaction(param_p2_r,param_p2_i,param_p2_p) + + r.add_reaction([("y",1)],[("h",1)],[("y",1)]) + r.add_reaction([("y",1)],[("Q",1)],[("yp",1)]) + r.add_reaction([("y",1),("yp",1)],[("h",1)],[("Y",1)]) + r.add_reaction([("z",1)],[("h",1)],[("z",1)]) + r.add_reaction([("z",1)], [("X",1)], [("zp",1)]) + r.add_reaction([("z",1),("zp",1)],[("h",1)],[("Z",1)]) + r.add_reaction([("U",1),("X",1)],[("h",1)],[("Q",1)]) + + c = ContextAutomatonWithConcentrations(r) + c.add_init_state("0") + c.add_state("1") + + # the experiments starts with adding x and y: + c.add_transition("0", [("x", 1),("y", 1)], "1") + + # for all the remaining steps we have empty context sequences + c.add_transition("1", [], "1") + + rc = ReactionSystemWithAutomaton(r, c) + rc.show() + + f_x1 = ltl_G(True, ltl_Implies( + exact_state(["x"], all_entities), + ltl_X(bag_And(bag_Not("Y"), bag_Not("Z"), bag_Not("h")), exact_state(["x", "xp"], all_entities)))) + f_x2 = ltl_G(True, ltl_Implies( + exact_state(["x", "xp"], all_entities), + ltl_X(bag_Not("h"), exact_state("X", all_entities)))) + + reach_xp = ltl_F(True, "xp") + reach_X = ltl_F(True, "X") + + smt_rsc = SmtCheckerRSCParam(rc, optimise=cmd_args.optimise) + + smt_rsc.check_rsltl(formulae_list=[f_x1, f_x2, reach_xp, reach_X]) + + # smt_rsc.check_rsltl(formula=f_x1, print_witness=True) + + +def trivial_param(): + + r = ReactionSystemWithConcentrationsParam() + r.add_bg_set_entity(("x", 3)) + r.add_bg_set_entity(("c", 3)) + r.add_bg_set_entity(("final", 1)) + + param_p1 = r.get_param("P1") + + r.add_reaction(param_p1, [("c", 1)], [("final", 1)]) + # r.add_reaction([("x", 1)], [("c", 1)], [("final", 1)]) c = ContextAutomatonWithConcentrations(r) c.add_init_state("0") c.add_state("1") - c.add_transition("0", [("y",1),("m",1),("x",1)], "1") - c.add_transition("1", [("x",1)], "1") - c.add_transition("1", [("x",1),("h",1)], "1") - c.add_transition("1", [("x",2)], "1") - c.add_transition("1", [("x",2),("h",1)], "1") - c.add_transition("1", [("h",1)], "1") + c.add_transition("0", [("x", 1)], "1") + c.add_transition("1", [], "1") + rc = ReactionSystemWithAutomaton(r, c) + rc.show() + smt_rsc = SmtCheckerRSCParam(rc) + + f1 = ltl_F(True, bag_entity("final") >= 1) + + # f1 = Formula_rsLTL.f_F( + # BagDescription.f_TRUE(), + # BagDescription.f_entity("final") >= 1) + + # + # WARNING: depth limit is set + # + smt_rsc.check_rsltl(formula=f1, max_level=10, print_witness=True) + + +def simple_param(cmd_args): + + r = ReactionSystemWithConcentrationsParam() + r.add_bg_set_entity(("x", 3)) + r.add_bg_set_entity(("y", 3)) + r.add_bg_set_entity(("c", 3)) + r.add_bg_set_entity(("z", 3)) + r.add_bg_set_entity(("final", 1)) + + param_p1 = r.get_param("P1") + + r.add_reaction([("x", 1)], [("c", 1)], [("y", 2)]) + # r.add_reaction([("y", 1)], [("c", 1)], [("z", 1)]) + r.add_reaction([("y", 1)], [("c", 1)], r.get_param("P2")) + r.add_reaction(param_p1, [("x", 1)], [("final", 1)]) + # r.add_reaction([("z", 1)], [("x", 1)], [("final", 1)]) + + c = ContextAutomatonWithConcentrations(r) + c.add_init_state("0") + c.add_state("1") + c.add_transition("0", [("x", 1)], "1") + c.add_transition("1", [], "1") + + rc = ReactionSystemWithAutomaton(r, c) + rc.show() + + smt_rsc = SmtCheckerRSCParam(rc, optimise=cmd_args.optimise) + + f1 = Formula_rsLTL.f_F( + BagDescription.f_TRUE(), + BagDescription.f_entity("final") >= 1) + + # + # WARNING: depth limit is set + # + smt_rsc.check_rsltl(formula=f1, max_level=10, print_witness=True, cont_if_sat=False) + + +def heat_shock_response_param(cmd_args, print_system=True): + + stress_temp = 42 + max_temp = 50 + + r = ReactionSystemWithConcentrationsParam() + r.add_bg_set_entity(("hsp",1)) + r.add_bg_set_entity(("hsf",1)) + r.add_bg_set_entity(("hsf2",1)) + r.add_bg_set_entity(("hsf3",1)) + r.add_bg_set_entity(("hse",1)) + r.add_bg_set_entity(("mfp",1)) + r.add_bg_set_entity(("prot",1)) + r.add_bg_set_entity(("hsf3:hse",1)) + r.add_bg_set_entity(("hsp:mfp",1)) + r.add_bg_set_entity(("hsp:hsf",1)) + r.add_bg_set_entity(("stress",1)) + r.add_bg_set_entity(("no_stress",1)) + + param_p1 = r.get_param("P1") + + r.add_reaction([("hsf",1)], [("hsp",1)], [("hsf3",1)]) + r.add_reaction([("hsf",1),("hsp",1),("mfp",1)], [], [("hsf3",1)]) + r.add_reaction([("hsf3",1)], [("hse",1),("hsp",1)],[("hsf",1)]) + r.add_reaction([("hsf3",1),("hsp",1),("mfp",1)], [("hse",1)], [("hsf",1)]) + r.add_reaction([("hsf3",1),("hse",1)], [("hsp",1)], [("hsf3:hse",1)]) + r.add_reaction([("hsp",1),("hsf3",1),("mfp",1),("hse",1)],[], [("hsf3:hse",1)]) + r.add_reaction([("hse",1)], [("hsf3",1)], [("hse",1)]) + r.add_reaction([("hsp",1),("hsf3",1),("hse",1)], [("mfp",1)], [("hse",1)]) + r.add_reaction([("hsf3:hse",1)], [("hsp",1)], [("hsp",1),("hsf3:hse",1)]) + r.add_reaction([("hsp",1),("mfp",1),("hsf3:hse",1)],[], [("hsp",1),("hsf3:hse",1)]) + r.add_reaction([("hsf",1),("hsp",1)], [("mfp",1)], [("hsp:hsf",1)]) + r.add_reaction([("hsp:hsf",1),("stress",1)],[], [("hsf",1),("hsp",1)]) + r.add_reaction([("hsp:hsf",1)], [("stress",1)],[("hsp:hsf",1)]) + r.add_reaction([("hsp",1),("hsf3:hse",1)], [("mfp",1)], [("hse",1),("hsp:hsf",1)]) + r.add_reaction([("stress",1),("prot",1)], [], [("mfp",1),("prot",1)]) + r.add_reaction([("prot",1)], [("stress",1)], [("prot",1)]) + r.add_reaction([("hsp",1),("mfp",1)], [], [("hsp:mfp",1)]) + r.add_reaction([("mfp",1)], [("hsp",1)], [("mfp",1)]) + r.add_reaction([("hsp:mfp",1)], [], [("hsp",1),("prot",1)]) + + # r.add_reaction_inc("temp", "heat", [("temp",1)],[("temp",max_temp)]) + # r.add_reaction_dec("temp", "cool", [("temp",1)],[]) + + # r.add_permanency("temp",[("heat",1),("cool",1)]) + + c = ContextAutomatonWithConcentrations(r) + c.add_init_state("0") + c.add_state("1") + c.add_transition("0", [("hsf",1),("prot",1),("hse",1),("no_stress",1)], "1") + c.add_transition("1", [("stress",1)], "1") + c.add_transition("1", [("no_stress",1)], "1") + c.add_transition("1", [], "1") rc = ReactionSystemWithAutomaton(r,c) + + if print_system: + rc.show() + + # prop_req = [("hsp:hsf",1),("hse",1),("prot",1)] + # prop_block = [("temp",stress_temp)] + # prop_req = [ ("mfp",1) ] + # prop_block = [ ] + # prop = (prop_req,prop_block) + # rs_prop = (state_translate_rsc2rs(prop_req),state_translate_rsc2rs(prop_block)) + # + + # f_reach_mfp = Formula_rsLTL.f_F(BagDescription.f_TRUE(), (BagDescription.f_entity("mfp") > 0) ) + f_reach_mfp = ltl_F(True, bag_entity("mfp") > 0) + + f_reach_hspmfp = Formula_rsLTL.f_F(BagDescription.f_TRUE(), (BagDescription.f_entity("hsp:mfp") > 0) ) + + smt_rsc = SmtCheckerRSCParam(rc, optimise=cmd_args.optimise) + + smt_rsc.check_rsltl(formulae_list=[f_reach_mfp, f_reach_hspmfp]) + + + # (1) if we keep increasing the temperature the protein will eventually misfold + # f_mfp_when_heating = Formula_rsLTL.f_X(BagDescription.f_TRUE(), + # Formula_rsLTL.f_F(BagDescription.f_entity("heat") > 0, (BagDescription.f_entity("mfp") > 0) ) + # ) + # smt_rsc.check_rsltl(formula=f_mfp_when_heating) + + # (2) when heating, we finally exceed the stress_temp + # f_2 = Formula_rsLTL.f_X(BagDescription.f_TRUE(), + # Formula_rsLTL.f_F(BagDescription.f_entity("heat") > 0, (BagDescription.f_entity("temp") > stress_temp)) + # ) + # smt_rsc.check_rsltl(formula=f_2) + + # smt_rsc.check_reachability(prop,max_level=40) + +def example44_param(): + + r = ReactionSystemWithConcentrationsParam() + r.add_bg_set_entity(("x", 2)) + r.add_bg_set_entity(("y", 4)) + r.add_bg_set_entity(("h", 2)) + r.add_bg_set_entity(("m", 1)) + + r.add_reaction([("y", 1), ("x", 1)], [("y", 2), ("h", 1)], [("y", 2)]) + r.add_reaction([("y", 2), ("x", 2)], [("y", 3), ("h", 1)], [("y", 3)]) + r.add_reaction([("y", 3), ("h", 1)], [("y", 4), ("h", 2)], [("y", 4)]) + r.add_reaction([("y", 4), ("h", 1)], [("h", 2)], [("y", 3)]) + r.add_reaction([("y", 4), ("x", 2)], [("h", 1)], [("y", 2)]) + r.add_reaction([("m", 1)], [("y", 3)], [("m", 1)]) + + c = ContextAutomatonWithConcentrations(r) + c.add_init_state("0") + c.add_state("1") + c.add_transition("0", [("y", 1), ("m", 1), ("x", 1)], "1") + c.add_transition("1", [("x", 1)], "1") + c.add_transition("1", [("x", 1), ("h", 1)], "1") + c.add_transition("1", [("x", 2)], "1") + c.add_transition("1", [("x", 2), ("h", 1)], "1") + c.add_transition("1", [("h", 1)], "1") + + rc = ReactionSystemWithAutomaton(r, c) rc.show() - smt_rsc = SmtCheckerRSC(rc) - + smt_rsc = SmtCheckerRSCParam(rc) + # Universal property which seems to be true: (holds also existentially) - f1 = Formula_rsLTL.f_G(BagDescription.f_entity("x") > 0, - Formula_rsLTL.f_Implies( - (BagDescription.f_entity('y') == 2), - Formula_rsLTL.f_X( - (BagDescription.f_entity("x") > 1), - BagDescription.f_entity("y") >= 3 - ) + f1 = Formula_rsLTL.f_G(BagDescription.f_entity("x") > 0, + Formula_rsLTL.f_Implies( + (BagDescription.f_entity('y') == 2), + Formula_rsLTL.f_X( + (BagDescription.f_entity("x") > 1), + BagDescription.f_entity("y") >= 3 ) ) - + ) + # lets see if we can find a counterexample to this property: - neg_f1 = Formula_rsLTL.f_F(BagDescription.f_entity("x") > 0, - Formula_rsLTL.f_And( - (BagDescription.f_entity('y') == 2), - Formula_rsLTL.f_X( - (BagDescription.f_entity("x") > 1), - BagDescription.f_entity("y") < 3 - ) + neg_f1 = Formula_rsLTL.f_F(BagDescription.f_entity("x") > 0, + Formula_rsLTL.f_And( + (BagDescription.f_entity('y') == 2), + Formula_rsLTL.f_X( + (BagDescription.f_entity("x") > 1), + BagDescription.f_entity("y") < 3 ) ) - + ) + # we fix the property f1 # this one holds: - f2 = Formula_rsLTL.f_G(BagDescription.f_entity("x") > 0, - Formula_rsLTL.f_Implies( - (BagDescription.f_entity('y') == 2), + f2 = Formula_rsLTL.f_G( + BagDescription.f_entity("x") > 0, Formula_rsLTL.f_Implies( + (BagDescription.f_entity('y') == 2), Formula_rsLTL.f_X( - (BagDescription.f_entity("x") > 1) & (BagDescription.f_entity("h") < 1), - BagDescription.f_entity("y") >= 3 - ) - ) - ) - + (BagDescription.f_entity("x") > 1) & + (BagDescription.f_entity("h") < 1), + BagDescription.f_entity("y") >= 3))) + # neg_f1 = Formula_rsLTL.f_X(BagDescription.f_TRUE(), Formula_rsLTL.f_F(BagDescription.f_entity("x") > 0, # Formula_rsLTL.f_And( # (BagDescription.f_entity('y') > 0), @@ -90,6 +331,83 @@ def example44(): # ) # )) smt_rsc.check_rsltl(formula=neg_f1) + + +def example44(): + + r = ReactionSystemWithConcentrations() + r.add_bg_set_entity(("x", 2)) + r.add_bg_set_entity(("y", 4)) + r.add_bg_set_entity(("h", 2)) + r.add_bg_set_entity(("m", 1)) + + r.add_reaction([("y", 1), ("x", 1)], [("y", 2), ("h", 1)], [("y", 2)]) + r.add_reaction([("y", 2), ("x", 2)], [("y", 3), ("h", 1)], [("y", 3)]) + r.add_reaction([("y", 3), ("h", 1)], [("y", 4), ("h", 2)], [("y", 4)]) + r.add_reaction([("y", 4), ("h", 1)], [("h", 2)], [("y", 3)]) + r.add_reaction([("y", 4), ("x", 2)], [("h", 1)], [("y", 2)]) + r.add_reaction([("m", 1)], [("y", 3)], [("m", 1)]) + + c = ContextAutomatonWithConcentrations(r) + c.add_init_state("0") + c.add_state("1") + c.add_transition("0", [("y", 1), ("m", 1), ("x", 1)], "1") + c.add_transition("1", [("x", 1)], "1") + c.add_transition("1", [("x", 1), ("h", 1)], "1") + c.add_transition("1", [("x", 2)], "1") + c.add_transition("1", [("x", 2), ("h", 1)], "1") + c.add_transition("1", [("h", 1)], "1") + + rc = ReactionSystemWithAutomaton(r, c) + rc.show() + smt_rsc = SmtCheckerRSC(rc) + + # Universal property which seems to be true: (holds also existentially) + f1 = Formula_rsLTL.f_G(BagDescription.f_entity("x") > 0, + Formula_rsLTL.f_Implies( + (BagDescription.f_entity('y') == 2), + Formula_rsLTL.f_X( + (BagDescription.f_entity("x") > 1), + BagDescription.f_entity("y") >= 3 + ) + ) + ) + + # lets see if we can find a counterexample to this property: + neg_f1 = Formula_rsLTL.f_F(BagDescription.f_entity("x") > 0, + Formula_rsLTL.f_And( + (BagDescription.f_entity('y') == 2), + Formula_rsLTL.f_X( + (BagDescription.f_entity("x") > 1), + BagDescription.f_entity("y") < 3 + ) + ) + ) + + # we fix the property f1 + # this one holds: + f2 = Formula_rsLTL.f_G( + BagDescription.f_entity("x") > 0, Formula_rsLTL.f_Implies( + (BagDescription.f_entity('y') == 2), + Formula_rsLTL.f_X( + (BagDescription.f_entity("x") > 1) & + (BagDescription.f_entity("h") < 1), + BagDescription.f_entity("y") >= 3))) + + # neg_f1 = Formula_rsLTL.f_X(BagDescription.f_TRUE(), Formula_rsLTL.f_F(BagDescription.f_entity("x") > 0, + # Formula_rsLTL.f_And( + # (BagDescription.f_entity('y') > 0), + # Formula_rsLTL.f_X( + # BagDescription.f_entity("x") > 0, + # Formula_rsLTL.f_X( + # BagDescription.f_entity("x") > 1, + # BagDescription.f_entity("y") < 3 + # ) + # ) + # ) + # )) + smt_rsc.check_rsltl(formula=neg_f1) + def heat_shock_response(print_system=True): diff --git a/rsltl_shortcuts.py b/rsltl_shortcuts.py new file mode 100644 index 0000000..1e3abcc --- /dev/null +++ b/rsltl_shortcuts.py @@ -0,0 +1,75 @@ +from logics import * + +##### SHORTCUTS + +def bag_True(): + return BagDescription.f_TRUE() + +def bag_entity(name): + return BagDescription.f_entity(name) + +def get_bag_if_str(arg): + if isinstance(arg, str): + return bag_entity(arg) > 0 + else: + return arg + +def bag_Not(a0): + a0 = get_bag_if_str(a0) + return BagDescription.f_Not(a0) + +def bag_And(*args): + assert len(args) > 1 + last = get_bag_if_str(args[0]) + for arg in args[1:]: + last = BagDescription.f_And( + last, get_bag_if_str(arg)) + return last + +def exact_state(contained_entities, all_entities): + expr = [] + for ent_str in all_entities: + ent = bag_entity(ent_str) + if ent_str in contained_entities: + expr.append(ent > 0) + else: + expr.append(ent == 0) + + if len(expr) > 0: + + last = expr[0] + for e in expr[1:]: + last = BagDescription.f_And(last, e) + return last + + else: + assert False + +def ltl_F(ctx_arg, a0): + a0 = get_bag_if_str(a0) + return Formula_rsLTL.f_F(ctx_arg, a0) + +def ltl_G(ctx_arg, a0): + a0 = get_bag_if_str(a0) + return Formula_rsLTL.f_G(ctx_arg, a0) + +def ltl_X(ctx_arg, a0): + a0 = get_bag_if_str(a0) + return Formula_rsLTL.f_X(ctx_arg, a0) + +def ltl_And(*args): + assert len(args) > 1 + last = get_bag_if_str(args[0]) + for arg in args[1:]: + last = Formula_rsLTL.f_And(last, get_bag_if_str(arg)) + return last + +def ltl_Not(a0): + a0 = get_bag_if_str(a0) + return Formula_rsLTL.f_Not(a0) + +def ltl_Implies(a0, a1): + a0 = get_bag_if_str(a0) + a1 = get_bag_if_str(a1) + return Formula_rsLTL.f_Implies(a0, a1) + diff --git a/rssmt.py b/rssmt.py index e37f75e..6e50d46 100755 --- a/rssmt.py +++ b/rssmt.py @@ -6,6 +6,8 @@ """ +import argparse + from rs import * from smt import * import sys @@ -14,38 +16,50 @@ import rs_testing from colour import * -import resource - profiling = False -################################################################## +if profiling: + import resource + +################################################################## version = "2017/04/08/00" rsmc_banner = """ Reaction Systems SMT-Based Model Checking Version: """ + version + """ -Author: Artur Męski / +Author: Artur Meski / """ -################################################################## +################################################################## + def print_banner(): print() for line in rsmc_banner.split("\n"): - print(colour_str(C_GREEN, " " + 3*"-" + " "), line) + print(colour_str(C_GREEN, " " + 3 * "-" + " "), line) print() -################################################################## +################################################################## + def main(): """Main function""" - - print_banner() - rs_testing.run_tests() -################################################################## - + parser = argparse.ArgumentParser() + + parser.add_argument("-v", "--verbose", + help="turn verbosity on", action="store_true") + parser.add_argument("-o", "--optimise", + help="minimise the parametric computation result", action="store_true") + + args = parser.parse_args() + + print_banner() + rs_testing.run_tests(args) + +################################################################## + if __name__ == "__main__": try: if profiling: diff --git a/smt/__init__.py b/smt/__init__.py index d4c1b71..8958447 100644 --- a/smt/__init__.py +++ b/smt/__init__.py @@ -1,5 +1,6 @@ #from smt.smt_checker import SmtChecker from smt.smt_checker_rs import SmtCheckerRS from smt.smt_checker_rsc import SmtCheckerRSC +from smt.smt_checker_rsc_param import SmtCheckerRSCParam from smt.smt_checker_rs_na import SmtCheckerRSNA from smt.smt_checker_distrib_rs import SmtCheckerDistribRS diff --git a/smt/smt_checker_rsc_param.py b/smt/smt_checker_rsc_param.py new file mode 100644 index 0000000..fb264ed --- /dev/null +++ b/smt/smt_checker_rsc_param.py @@ -0,0 +1,961 @@ +""" +SMT-based Model Checking Module for RS with Concentrations and Context Automaton +""" + +from z3 import * +from time import time +from sys import stdout +from itertools import chain +import resource +from colour import * + +from logics import rsLTL_Encoder + +from rs.reaction_system_with_concentrations_param import ParameterObj, is_param + +# def simplify(x): +# return x + +def z3_max(a, b): + return If(a > b, a, b) + +class SmtCheckerRSCParam(object): + + def __init__(self, rsca, optimise=False): + + rsca.sanity_check() + + if not rsca.is_concentr_and_param_compatible(): + raise RuntimeError("RS and CA with concentrations (and parameters) expected") + + self.rs = rsca.rs + self.ca = rsca.ca + + self.optimise = optimise + + self.initialise() + + def initialise(self): + """Initialises all the variables used by the checker""" + + ### "Currently" used variables (loaded from self.path_v...) + self.v = None + self.v_ctx = None + self.ca_state = None + + # intermediate products: + self.v_improd = None + self.v_improd_for_entities = None + + ### Per-path variables + self.path_v = dict() + self.path_v_ctx = dict() + self.path_ca_state = dict() + + # intermediate products: + self.path_v_improd = dict() + self.path_v_improd_for_entities = dict() + + # parameters: + self.v_param = dict() + + self.next_level_to_encode = 0 + + # this is probably not needed anymore: + # self.producible_entities = self.rs.get_producible_entities() + # self.improducible_entities = set(self.rs.get_state_ids(self.rs.background_set)) - self.producible_entities + + # WARNING: improd vs. improducible + # there is some confusion related to the variables naming: + # improd - intermediate products + # improducible - entities that are never produces (there is no reaction that produces that entity) + + # TODO: number of loops == number of paths + self.loop_position = None + self.path_loop_position = dict() + + if self.optimise: + self.solver = Optimize() + else: + self.solver = Solver() #For("QF_FD") + + self.verification_time = None + + self.prepare_param_variables() + + def reset(self): + """Reinitialises the state of the checker""" + + self.initialise() + + def prepare_all_variables(self, num_of_paths): + + for path_idx in range(num_of_paths): + self.prepare_all_path_variables(path_idx) + self.next_level_to_encode += 1 + + def prepare_all_path_variables(self, path_idx): + """Prepares the variables for a given path index""" + + print_info("Preparing variables for path={:d} (level={:d})".format(path_idx, self.next_level_to_encode)) + + self.prepare_state_variables(path_idx) + self.prepare_context_variables(path_idx) + self.prepare_intermediate_product_variables(path_idx) + self.prepare_loop_position_variables(path_idx) + + def prepare_loop_position_variables(self, path_idx): + """Prepares the variables for loop positions""" + + self.path_loop_position[path_idx] = Int("p{:d}_loop_pos".format(path_idx)) + + def prepare_context_variables(self, path_idx): + """Prepares all the context variables""" + + level = self.next_level_to_encode + + self.path_v_ctx.setdefault(path_idx, []) + assert len(self.path_v_ctx[path_idx]) == level + + variables = [] + for entity in self.rs.background_set: + new_var = Int("p{:d}C{:d}_{:s}".format(path_idx, level, entity)) + variables.append(new_var) + + self.path_v_ctx[path_idx].append(variables) + + def prepare_state_variables(self, path_idx): + """Prepares all the state variables""" + + level = self.next_level_to_encode + + # RS vars + self.path_v.setdefault(path_idx, []) + assert len(self.path_v[path_idx]) == level + + variables = [] + for entity in self.rs.background_set: + new_var = Int("p{:d}L{:d}_{:s}".format(path_idx, level, entity)) + variables.append(new_var) + self.path_v[path_idx].append(variables) + + # Context automaton states: + self.path_ca_state.setdefault(path_idx, []) + assert len(self.path_ca_state[path_idx]) == level + + ca_state_var = Int("p{:d}CA{:d}_state".format(path_idx, level)) + self.path_ca_state[path_idx].append(ca_state_var) + + def prepare_intermediate_product_variables(self, path_idx): + """ + Prepares the intermediate product variables + carrying the individual concentration levels produced + by the reactions. + + These variables are used later on to encode the final + concentration levels for all the entities + """ + + level = self.next_level_to_encode + + self.path_v_improd.setdefault(path_idx, []) + self.path_v_improd_for_entities.setdefault(path_idx, []) + + assert len(self.path_v_improd[path_idx]) == level + assert len(self.path_v_improd_for_entities[path_idx]) == level + + if level < 1: + # + # If we are at level==0, we add a dummy "level" + # to match the indices of of the successors + # which are always at level+1. + # + self.path_v_improd[path_idx].append(None) + self.path_v_improd_for_entities[path_idx].append(None) + else: + + reactions_dict = dict() + number_of_reactions = len(self.rs.reactions) + + all_entities_dict = dict() + + for reaction in self.rs.reactions: + *_, products = reaction + + reaction_id = self.rs.reactions.index(reaction) + + entities_dict = dict() + + if is_param(products): + + for entity in self.rs.set_of_bgset_ids: + entity_name = self.rs.get_entity_name(entity) + new_var = Int("p{:d}L{:d}_ImProd_r{:d}_{:s}".format(path_idx, level, reaction_id, entity_name)) + entities_dict[entity] = new_var + + all_entities_dict.setdefault(entity, []) + all_entities_dict[entity].append(new_var) + + else: + + for entity, conc in products: + entity_name = self.rs.get_entity_name(entity) + new_var = Int("p{:d}L{:d}_ImProd_r{:d}_{:s}".format(path_idx, level, reaction_id, entity_name)) + entities_dict[entity] = new_var + + all_entities_dict.setdefault(entity, []) + all_entities_dict[entity].append(new_var) + + reactions_dict[reaction_id] = entities_dict + + self.path_v_improd[path_idx].append(reactions_dict) + self.path_v_improd_for_entities[path_idx].append(all_entities_dict) + + def prepare_param_variables(self): + """ + Prepares variables for parameters + + A parameter (it's valuation) is a subset of the background set, + therefore we need separate variables for each element of the + background set. + """ + + for param_name in self.rs.parameters.keys(): + + # we start collecting bg-related vars for the given param + vars_for_param = [] + + for entity in self.rs.set_of_bgset_ids: + new_var = Int("Pm{:s}_{:s}".format(param_name, self.rs.get_entity_name(entity))) + vars_for_param.append(new_var) + + self.v_param[param_name] = vars_for_param + + def load_varset_for_path(self, path_idx): + """ + Loads the the variables for the path with path_idx + """ + self.v = self.path_v[path_idx] + self.v_ctx = self.path_v_ctx[path_idx] + self.ca_state = self.path_ca_state[path_idx] + self.v_improd = self.path_v_improd[path_idx] + self.v_improd_for_entities = self.path_v_improd_for_entities[path_idx] + self.loop_position = self.path_loop_position[path_idx] + + def enc_param_concentration_levels_assertion(self): + """ + Assertions for the parameter variables + """ + + if len(self.v_param) == 0: + return True + + enc_param_gz = True + enc_non_empty = True + + for param_vars in self.v_param.values(): + + enc_param_at_least_one = False + for pvar in param_vars: + + # TODO: fixed upper limit: 100 (have a per-param setting for that) + enc_param_gz = simplify(And(enc_param_gz, pvar >= 0, pvar < 100)) + enc_param_at_least_one = simplify(Or(enc_param_at_least_one, pvar > 0)) + + enc_non_empty = simplify(And(enc_non_empty, enc_param_at_least_one)) + + return simplify(And(enc_param_gz, enc_non_empty)) + + def assert_param_optimisation(self): + + for param_vars in self.v_param.values(): + for pvar in param_vars: + self.solver.add_soft(pvar < 1) + + def enc_concentration_levels_assertion(self, level, path_idx): + """ + Encodes assertions that (some) variables need to be >=0 + + We do not need to actually control all the variables, + only those that can possibly go below 0. + """ + + print_info("Concentration level assertions for path={:d} (level={:d})".format(path_idx, level)) + + enc_gz = True + + for e_i in self.rs.set_of_bgset_ids: + var = self.path_v[path_idx][level][e_i] + var_ctx = self.path_v_ctx[path_idx][level][e_i] + e_max = self.rs.get_max_concentration_level(e_i) + enc_gz = simplify(And(enc_gz, var >= 0, var_ctx >= 0, var <= e_max, var_ctx <= e_max)) + + vars_per_reaction = self.path_v_improd_for_entities[path_idx][level + 1] + if e_i in vars_per_reaction: + for var_improd in vars_per_reaction[e_i]: + enc_gz = simplify(And(enc_gz, var_improd >= 0, var_improd <= e_max)) + + return enc_gz + + def enc_init_state(self, level, path_idx): + """Encodes the initial state at the given level""" + + rs_init_state_enc = True + + for v in self.path_v[path_idx][level]: + # the initial concentration levels are zeroed + rs_init_state_enc = simplify(And(rs_init_state_enc, v == 0)) + + ca_init_state_enc = self.path_ca_state[path_idx][level] == self.ca.get_init_state_id() + + init_state_enc = simplify(And(rs_init_state_enc, ca_init_state_enc)) + + return init_state_enc + + def enc_transition_relation(self, level, path_idx): + return simplify( + And(self.enc_rs_trans(level, path_idx), + self.enc_automaton_trans(level, path_idx))) + + def enc_single_reaction(self, level, path_idx, reaction): + """ + Encodes a single reaction + + For encoding the products we use intermediate variables: + + * each reaction has its own product variables, + + * those are meant to be used to compute the MAX concentration + + """ + + reactants, inhibitors, products = reaction + + # we need reaction_id to find the intermediate product variable + reaction_id = self.rs.reactions.index(reaction) + + # ** REACTANTS ******************************************* + enc_reactants = True + if is_param(reactants): + param_name = reactants.name + for entity in self.rs.set_of_bgset_ids: + enc_reactants = And(enc_reactants, Or( + self.v_param[param_name][entity] == 0, + self.path_v[path_idx][level][entity] >= self.v_param[param_name][entity], + self.path_v_ctx[path_idx][level][entity] >= self.v_param[param_name][entity])) + else: + for entity, conc in reactants: + enc_reactants = And(enc_reactants, Or( + self.path_v[path_idx][level][entity] >= conc, + self.path_v_ctx[path_idx][level][entity] >= conc)) + + # ** INHIBITORS ****************************************** + enc_inhibitors = True + if is_param(inhibitors): + param_name = inhibitors.name + for entity in self.rs.set_of_bgset_ids: + enc_inhibitors = And(enc_inhibitors, Or( + self.v_param[param_name][entity] == 0, + And( + self.path_v[path_idx][level][entity] < self.v_param[param_name][entity], + self.path_v_ctx[path_idx][level][entity] < self.v_param[param_name][entity]))) + else: + for entity, conc in inhibitors: + enc_inhibitors = And(enc_inhibitors, And( + self.path_v[path_idx][level][entity] < conc, + self.path_v_ctx[path_idx][level][entity] < conc)) + + # ** PRODUCTS ******************************************* + enc_products = True + if is_param(products): + param_name = products.name + for entity in self.rs.set_of_bgset_ids: + enc_products = simplify(And(enc_products, + self.path_v_improd[path_idx][level + 1][reaction_id][entity] == self.v_param[param_name][entity])) + else: + for entity, conc in products: + enc_products = simplify(And(enc_products, + self.path_v_improd[path_idx][level + 1][reaction_id][entity] == conc)) + + # Nothing is produced (when the reaction is disabled) + enc_no_prod = True + if is_param(products): + for entity in self.rs.set_of_bgset_ids: + enc_no_prod = And(enc_no_prod, + self.path_v_improd[path_idx][level + 1][reaction_id][entity] == 0) + else: + for entity, _ in products: + enc_no_prod = simplify(And(enc_no_prod, + self.path_v_improd[path_idx][level + 1][reaction_id][entity] == 0)) + + # + # (R and I) iff P + # + enc_enabled = And(enc_reactants, enc_inhibitors) == enc_products + + # + # ~(R and I) iff P_zero + # + enc_not_enabled = Not(And(enc_reactants, enc_inhibitors)) == enc_no_prod + + enc_reaction = And(enc_enabled, enc_not_enabled) + + return enc_reaction + + def enc_general_reaction_enabledness(self, level, path_idx): + """ + General enabledness condition for reactions + + The necessary condition for a reaction to be enabled is + that the state is not empty, i.e., at least one entity + is present in the current state. + + This condition must be used when there are parametric reactions + because parameters could have all the entities set to zero and + that immediately allows for all the conditions on the reactants + to be fulfilled: (entity <= param) -> (0 <= 0) + """ + + enc_cond = False + for entity in self.rs.set_of_bgset_ids: + enc_cond = simplify(Or(enc_cond, self.path_v[path_idx][level][entity] > 0, self.path_v_ctx[path_idx][level][entity] > 0)) + + return enc_cond + + def enc_rs_trans(self, level, path_idx): + """Encodes the transition relation""" + + # + # IMPORTANT NOTE + # + # We need to make sure we do something about the UNUSED ENTITIES + # that is, those that are never produced. + # + # They should have concentration levels set to 0. + # + # That needs to happen automatically (in the MAX encoding) -- for the parametric + # case it makes no sense to identify the entities that are never produced, unless + # we have no parameters as products (special case, so that could be an + # optimisation) + # + + enc_trans = True + + for reaction in self.rs.reactions: + enc_reaction = self.enc_single_reaction(level, path_idx, reaction) + enc_trans = simplify(And(enc_trans, enc_reaction)) + + # Next we encode the MAX concentration values: + # we collect those from the intermediate product variables + + enc_max_prod = True + + # Save all the intermediate product variables for a given level: + # + # - Intermediate products of (level+1) correspond to the next level + # + # {reactants & inhibitors}[level] + # => + # {improd}[level+1] + # => + # {products}[level+1] + # + current_v_improd_for_entities = self.path_v_improd_for_entities[path_idx][level + 1] + for entity in self.rs.set_of_bgset_ids: + per_reaction_vars = current_v_improd_for_entities.get(entity, []) + enc_max_prod = simplify( + And(enc_max_prod, self.path_v[path_idx][level + 1][entity] == self.enc_max(per_reaction_vars))) + + # make sure at least one entity is >0 + enc_general_cond = self.enc_general_reaction_enabledness(level, path_idx) + + enc_trans_with_max = simplify(And(enc_general_cond, enc_max_prod, enc_trans)) + + # print(enc_trans_with_max) + + return enc_trans_with_max + + def enc_max(self, elements): + + enc = None + + if elements == []: + enc = 0 + + elif len(elements) == 1: + enc = z3_max(0, elements[0]) + + elif len(elements) > 1: + + enc = 0 + for i in range(len(elements) - 1): + enc = z3_max(enc, z3_max(elements[i], elements[i + 1])) + + return enc + + def enc_automaton_trans(self, level, path_idx): + """Encodes the transition relation for the context automaton""" + + enc_trans = False + + for src,ctx,dst in self.ca.transitions: + src_enc = self.path_ca_state[path_idx][level] == src + dst_enc = self.path_ca_state[path_idx][level+1] == dst + + all_ent = set(range(len(self.rs.background_set))) + + incl_ctx = set([e for e,c in ctx]) + excl_ctx = all_ent - incl_ctx + + ctx_enc = True + + for e,c in ctx: + ctx_enc = simplify(And(ctx_enc, self.path_v_ctx[path_idx][level][e] == c)) + + for e in excl_ctx: + ctx_enc = simplify(And(ctx_enc, self.path_v_ctx[path_idx][level][e] == 0)) + + cur_trans = simplify(And(src_enc, ctx_enc, dst_enc)) + enc_trans = simplify(Or(enc_trans, cur_trans)) + + return enc_trans + + def enc_exact_state(self, level, state): + """Encodes the state at the given level with the exact concentration values""" + + raise RuntimeError("Should not be used with RSC") + + def enc_min_state(self, level, state): + """Encodes the state at the given level with the minimal required concentration levels""" + + enc = True + for ent,conc in state: + e_id = self.rs.get_entity_id(ent) + enc = And(enc, self.v[level][e_id] >= conc) + + # state_ids = self.rs.get_state_ids(state) + # + # for entity in state_ids: + # enc = And(enc, self.v[level][entity]) + + return simplify(enc) + + def enc_state_with_blocking(self, level, prop): + """Encodes the state at the given level with blocking certain concentrations""" + + required,blocked = prop + + enc = True + for ent,conc in required: + e_id = self.rs.get_entity_id(ent) + enc = And(enc, self.v[level][e_id] >= conc) + + for ent,conc in blocked: + e_id = self.rs.get_entity_id(ent) + enc = And(enc, self.v[level][e_id] < conc) + + return simplify(enc) + + def decode_witness(self, max_level, path_idx, print_model=False): + """ + Decodes the witness + + Also decodes the parameters + """ + + m = self.solver.model() + + if print_model: + print(m) + + for level in range(max_level + 1): + + print("\n{: >70}".format("[ level=" + repr(level) + " ]")) + + print(" State: {", end=""), + for var_id in range(len(self.path_v[path_idx][level])): + var_rep = repr(m[self.path_v[path_idx][level][var_id]]) + if not var_rep.isdigit(): + raise RuntimeError( + "unexpected: representation is not a positive integer") + if int(var_rep) > 0: + print( + " " + self.rs.get_entity_name(var_id) + "=" + var_rep, + end="") + # print(" " + repr(m[self.v[level][var_id]]), end="") + print(" }") + + if level != max_level: + print(" Context set: ", end="") + print("{", end="") + for var_id in range(len(self.path_v[path_idx][level])): + var_rep = repr(m[self.path_v_ctx[path_idx][level][var_id]]) + if not var_rep.isdigit(): + raise RuntimeError( + "unexpected: representation is not a positive integer") + if int(var_rep) > 0: + print( + " " + self.rs.get_entity_name(var_id) + "=" + var_rep, + end="") + print(" }") + + # Parameters + + print("\n\n Parameters:\n") + for param_name in self.rs.parameters.keys(): + print("{: >6}: ".format(param_name), end="") + print("{", end="") + + params = self.v_param[param_name] + + for entity in self.rs.set_of_bgset_ids: + var_rep = repr(m[params[entity]]) + if not var_rep.isdigit(): + raise RuntimeError( + "unexpected: representation is not a positive integer") + if int(var_rep) > 0: + print( + " " + str(self.rs.get_entity_name(entity)) + "=" + str(var_rep), + end="") + print(" }") + + print("\n") + + def check_rsltl( + self, formulae_list, + print_witness=True, + print_time=True, print_mem=True, + max_level=None, cont_if_sat=False): + """ + Bounded Model Checking for rsLTL properties + + * print_witness -- prints the decoded witness + * print_time -- prints the time consumed + * print_mem -- prints the memory consumed + * max_level -- if not None, the methods + stops at the specified level + * cont_if_sat -- if True, then the method + continues up until max_level is + reached (even if sat found) + """ + + if not isinstance(formulae_list, (list, tuple)): + print_error("Expected a list of formulae") + + self.reset() + + num_of_paths = len(formulae_list) + + # formula = formulae[0] + + print_info("Running rsLTL bounded model checking") + print_info("Tested formulae:") + for form in formulae_list: + print_info(" "*4 + str(form)) + + if print_time: + start = resource.getrusage(resource.RUSAGE_SELF).ru_utime + + self.prepare_all_variables(num_of_paths) + + self.load_varset_for_path(0) + + # initial states for all the paths + initial_states = [] + for path_idx in range(num_of_paths): + initial_states.append(self.enc_init_state(0, path_idx)) + self.solver_add(initial_states) + + self.current_level = 0 + + self.prepare_all_variables(num_of_paths) + + # assertions for all the paths and parameters + additional_assertions = [] + for path_idx in range(num_of_paths): + additional_assertions.append(self.enc_concentration_levels_assertion(0, path_idx)) + additional_assertions.append(self.enc_param_concentration_levels_assertion()) + self.solver_add(additional_assertions) + + if self.optimise: + self.assert_param_optimisation() + + encoder = rsLTL_Encoder(self) + + print_info("Iterating...") + + while True: + + + + print( + "\n{:-^70}".format("[ Working at level=" + str(self.current_level) + " ]")) + # stdout.flush() + + # reachability test: + self.solver.push() + + enc_form = [] + for formula in formulae_list: + + path_idx = formulae_list.index(formula) + + print_info("Generating the encoding for {:s} ({:d} of {:d})".format( + str(formula), path_idx+1, len(formulae_list))) + + encoder.load_variables( + var_rs=self.path_v[path_idx], + var_ctx=self.path_v_ctx[path_idx], + var_loop_pos=self.path_loop_position[path_idx]) + + enc_form.append(encoder.get_encoding(formula, self.current_level)) + ncalls = encoder.get_ncalls() + + print_info("Cache hits: {:d}, encode calls: {:d} (approx: {:d})".format( + encoder.get_cache_hits(), ncalls[0], ncalls[1])) + + encoder.flush_cache() + + print_info("Adding the formulae to the solver...") + + # print (enc_form) + + self.solver_add(enc_form) + + print_info("Adding the loops encoding...") + self.solver_add(self.get_loop_encodings()) + + print_info("Testing satisfiability...") + result = self.solver.check() + if result == sat: + print_positive(green_str( + "SAT at level={:d}".format(self.current_level))) + # print(self.solver.model()) + if print_witness: + for formula in formulae_list: + path_idx = formulae_list.index(formula) + print("\n{:=^70}".format("[ WITNESS ]")) + print("\n Witness for: {:s}".format(str(formula))) + self.decode_witness(self.current_level, path_idx) + if not cont_if_sat: + break + else: + print_info("Unsat") + self.solver.pop() + + self.prepare_all_variables(num_of_paths) + + # assertions for all the paths + additional_assertions = [] + for path_idx in range(num_of_paths): + additional_assertions.append( + self.enc_concentration_levels_assertion(self.current_level + 1, path_idx)) + self.solver_add(additional_assertions) + + print_info("Unrolling the transition relation") + for path_idx in range(num_of_paths): + self.solver_add(self.enc_transition_relation(self.current_level, path_idx)) + + print( + "{:->70}".format("[ level=" + str(self.current_level) + " done ]")) + self.current_level += 1 + + if not max_level is None and self.current_level > max_level: + print("Stopping at level=" + str(max_level)) + break + + if print_time: + # stop = time() + stop = resource.getrusage(resource.RUSAGE_SELF).ru_utime + self.verification_time = stop - start + print() + print( + "\n[i] {: >60}".format( + " Time: " + repr(self.verification_time) + " s")) + + if print_mem: + print( + "[i] {: >60}".format( + " Memory: " + + repr( + resource.getrusage(resource.RUSAGE_SELF).ru_maxrss / + (1024 * 1024)) + " MB")) + + def dummy_unroll(self, levels): + """Unrolls the variables for testing purposes""" + + self.current_level = -1 + for i in range(levels+1): + self.prepare_all_variables() + self.current_level += 1 + + print(C_MARK_INFO + " Dummy Unrolling done.") + + def state_equality(self, level_A, level_B): + """Encodes equality of two states at two different levels""" + + eq_enc = True + + for e_i in range(len(self.rs.background_set)): + e_i_equality = self.v[level_A][e_i] == self.v[level_B][e_i] + eq_enc = simplify(And(eq_enc, e_i_equality)) + + eq_enc_ctxaut = self.ca_state[level_A] == self.ca_state[level_B] + eq_enc = simplify(And(eq_enc, eq_enc_ctxaut)) + + return eq_enc + + def get_loop_encodings(self): + + k = self.current_level + loop_var = self.loop_position + + loop_enc = True + + """ + (loop_var == i) means that there is a loop taking back to the state (i-1) + + Therefore, the encoding starts at 1, not at 0. + """ + + for i in range(1,k+1): + loop_enc = simplify(And(loop_enc, Implies( loop_var == i, self.state_equality(i-1, k) ))) + + return loop_enc + + def solver_add(self, expression): + """ + This is a solver.add() wrapper + """ + + if expression == True: + return + + if expression == False: + raise RuntimeError("Trying to assert False.") + + self.solver.add(expression) + + def check_reachability(self, state, print_witness=True, + print_time=True, print_mem=True, max_level=1000): + """Main testing function""" + + self.reset() + + if print_time: + # start = time() + start = resource.getrusage(resource.RUSAGE_SELF).ru_utime + + self.prepare_all_variables() + self.solver_add(self.enc_init_state(0)) + self.current_level = 0 + + self.prepare_all_variables() + + self.solver_add(self.enc_concentration_levels_assertion(0)) + + while True: + self.prepare_all_variables() + self.solver_add(self.enc_concentration_levels_assertion(self.current_level+1)) + + print("\n{:-^70}".format("[ Working at level=" + str(self.current_level) + " ]")) + stdout.flush() + + # reachability test: + print("[" + colour_str(C_BOLD, "i") + "] Adding the reachability test...") + self.solver.push() + + self.solver_add(self.enc_state_with_blocking(self.current_level,state)) + + result = self.solver.check() + if result == sat: + print("[" + colour_str(C_BOLD, "+") + "] " + colour_str(C_GREEN, "SAT at level=" + str(self.current_level))) + if print_witness: + print("\n{:=^70}".format("[ WITNESS ]")) + self.decode_witness(self.current_level) + break + else: + self.solver.pop() + + print("[" + colour_str(C_BOLD, "i") + "] Unrolling the transition relation") + self.solver_add(self.enc_transition_relation(self.current_level)) + + print("{:->70}".format("[ level=" + str(self.current_level) + " done ]")) + self.current_level += 1 + + if self.current_level > max_level: + print("Stopping at level=" + str(max_level)) + break + + if print_time: + # stop = time() + stop = resource.getrusage(resource.RUSAGE_SELF).ru_utime + self.verification_time = stop-start + print() + print("\n[i] {: >60}".format(" Time: " + repr(self.verification_time) + " s")) + + if print_mem: + print("[i] {: >60}".format(" Memory: " + repr(resource.getrusage(resource.RUSAGE_SELF).ru_maxrss/(1024*1024)) + " MB")) + + + def get_verification_time(self): + return self.verification_time + + + def show_encoding(self, state, print_witness=True, + print_time=False, print_mem=False, max_level=100): + """Encoding debug function""" + + self.reset() + + self.prepare_all_variables() + init_s = self.enc_init_state(0) + print(init_s) + self.solver_add(init_s) + self.current_level = 0 + + self.prepare_all_variables() + + while True: + self.prepare_all_variables() + + print("-----[ Working at level=" + str(self.current_level) + " ]-----") + stdout.flush() + + # reachability test: + print("[i] Adding the reachability test...") + self.solver.push() + + s = self.enc_min_state(self.current_level,state) + print("Test: ", s) + + self.solver_add(s) + + result = self.solver.check() + if result == sat: + print("\n[+] " + colour_str(C_RED, "SAT at level=" + str(self.current_level))) + if print_witness: + self.decode_witness(self.current_level) + break + else: + self.solver.pop() + + print("[i] Unrolling the transition relation") + t = self.enc_transition_relation(self.current_level) + print(t) + self.solver_add(t) + + print("-----[ level=" + str(self.current_level) + " done ]") + self.current_level += 1 + + if self.current_level > max_level: + print("Stopping at level=" + str(max_level)) + break + else: + x=input("Next level? ") + x=x.lower() + if not (x == "y" or x == "yes"): + break + + +# EOF \ No newline at end of file