From 55ff688fd03c872ac490c0de680fbea537eee761 Mon Sep 17 00:00:00 2001 From: Artur Meski Date: Wed, 9 Aug 2017 20:57:52 +0100 Subject: [PATCH 01/35] add RSC with param for the new encoding suitable for parametric verification --- rs/__init__.py | 1 + ...action_system_with_concentrations_param.py | 414 ++++++++++++++++++ 2 files changed, 415 insertions(+) create mode 100644 rs/reaction_system_with_concentrations_param.py diff --git a/rs/__init__.py b/rs/__init__.py index b500f08..01001bf 100644 --- a/rs/__init__.py +++ b/rs/__init__.py @@ -2,6 +2,7 @@ from rs.reaction_system import ReactionSystem from rs.context_automaton import ContextAutomaton from rs.reaction_system_with_concentrations import ReactionSystemWithConcentrations +from rs.reaction_system_with_concentrations_param import ReactionSystemWithConcentrationsParam from rs.context_automaton_with_concentrations import ContextAutomatonWithConcentrations from rs.extended_context_automaton import ExtendedContextAutomaton diff --git a/rs/reaction_system_with_concentrations_param.py b/rs/reaction_system_with_concentrations_param.py new file mode 100644 index 0000000..5c18e0f --- /dev/null +++ b/rs/reaction_system_with_concentrations_param.py @@ -0,0 +1,414 @@ +from sys import exit +from colour import * + +from rs.reaction_system import ReactionSystem + +class ReactionSystemWithConcentrationsParam(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:") + print(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 + 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): + + 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) + + +# 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() From e1db24c8a976e1b15f5de6f05bb13dba06548d9a Mon Sep 17 00:00:00 2001 From: Artur Meski Date: Sun, 13 Aug 2017 12:26:43 +0100 Subject: [PATCH 02/35] SMTChecker for RSC, for the new parametic-compatible encoding --- smt/smt_checker_rsc_param.py | 616 +++++++++++++++++++++++++++++++++++ 1 file changed, 616 insertions(+) create mode 100644 smt/smt_checker_rsc_param.py diff --git a/smt/smt_checker_rsc_param.py b/smt/smt_checker_rsc_param.py new file mode 100644 index 0000000..4604f87 --- /dev/null +++ b/smt/smt_checker_rsc_param.py @@ -0,0 +1,616 @@ +""" +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]: + rs_init_state_enc = simplify(And(rs_init_state_enc, v == 0)) # the initial concentration levels are zeroed + + 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 == []: + return simplify(self.v[level+1][prod_entity] == 0) # this should never happen + + 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)) + + # for reactants,inhibitors,products in rcts_for_prod_entity: + # enc_reactants = True + # enc_inhibitors = True + # # enc_products -- below + # + # for reactant,concentration in reactants: + # enc_reactants = simplify(And(enc_reactants, + # Or(self.v[level][reactant] >= concentration, self.v_ctx[level][reactant] >= concentration))) + # for inhibitor,concentration in inhibitors: + # enc_inhibitors = simplify(And(enc_inhibitors, + # And(self.v[level][inhibitor] < concentration, self.v_ctx[level][inhibitor] < concentration))) + # + # enc_products = self.v[level+1][products[0][0]] == products[0][1] + # + # enc_enabledness = simplify(Or(enc_enabledness, And(enc_reactants, enc_inhibitors))) + # enc_rct_prod = simplify(Or(enc_rct_prod, And(enc_reactants, enc_inhibitors, enc_products))) + + # -------- 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_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): + 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") + + # enc = True + # used_entities_ids = self.rs.get_state_ids(state) + # + # for ent,conc in state: + # e_id = self.rs.get_entity_id(ent) + # enc = And(enc, self.v[level][e_id] == conc) + # + # not_in_state = set(range(len(self.rs.background_set))) + # not_in_state = not_in_state.difference(set(used_entities_ids)) + # + # for entity in not_in_state: + # enc = And(enc, self.v[level][entity] == 0) + # + # return simplify(enc) + + 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, 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) + + 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 + From 27c86497e4d7096ee6f75f49b5e65fa64d6fdc1d Mon Sep 17 00:00:00 2001 From: Artur Meski Date: Sun, 13 Aug 2017 13:48:42 +0100 Subject: [PATCH 03/35] RSC with Param; RSCA for Param, SMTChecker for Param --- rs/reaction_system_with_automaton.py | 12 +++ ...action_system_with_concentrations_param.py | 85 ++++++++++++------- rs_testing.py | 81 +++++++++++++++++- smt/__init__.py | 1 + smt/smt_checker_rsc_param.py | 8 +- 5 files changed, 151 insertions(+), 36 deletions(-) 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_param.py b/rs/reaction_system_with_concentrations_param.py index 5c18e0f..20e2134 100644 --- a/rs/reaction_system_with_concentrations_param.py +++ b/rs/reaction_system_with_concentrations_param.py @@ -102,7 +102,10 @@ class ReactionSystemWithConcentrationsParam(ReactionSystem): return reactants,inhibitors,products def add_reaction(self, R, I, P): - """Adds a reaction""" + """Adds a reaction + + R, I, and P are sets of entities (not their IDs) + """ if P == []: raise RuntimeError("No products defined") @@ -201,6 +204,8 @@ class ReactionSystemWithConcentrationsParam(ReactionSystem): 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 @@ -247,6 +252,9 @@ class ReactionSystemWithConcentrationsParam(ReactionSystem): return reactions_by_prod def get_reaction_system(self): + """ + Translates RSC into RS + """ rs = ReactionSystem() @@ -372,36 +380,51 @@ class ReactionSystemWithConcentrationsParam(ReactionSystem): 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) - +# 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): # diff --git a/rs_testing.py b/rs_testing.py index 6ecb83c..4e61f7b 100644 --- a/rs_testing.py +++ b/rs_testing.py @@ -11,7 +11,86 @@ def run_tests(): # process() # heat_shock_response() # scalable_chain(print_system=True) - example44() + # example44() + example44_param() + +def example44_param(): + + r = ReactionSystemWithConcentrationsParam() + r.add_bg_set_entity(("x",2)) + r.add_bg_set_entity(("y",4)) + r.add_bg_set_entity(("h",2)) + r.add_bg_set_entity(("m",1)) + + r.add_reaction([("y",1),("x",1)], [("y",2),("h",1)], [("y",2)]) + r.add_reaction([("y",2),("x",2)], [("y",3),("h",1)], [("y",3)]) + r.add_reaction([("y",3),("h",1)], [("y",4),("h",2)], [("y",4)]) + r.add_reaction([("y",4),("h",1)], [("h",2)], [("y",3)]) + r.add_reaction([("y",4),("x",2)], [("h",1)], [("y",2)]) + r.add_reaction([("m",1)], [("y",3)], [("m",1)]) + + c = ContextAutomatonWithConcentrations(r) + c.add_init_state("0") + c.add_state("1") + c.add_transition("0", [("y",1),("m",1),("x",1)], "1") + c.add_transition("1", [("x",1)], "1") + c.add_transition("1", [("x",1),("h",1)], "1") + c.add_transition("1", [("x",2)], "1") + c.add_transition("1", [("x",2),("h",1)], "1") + c.add_transition("1", [("h",1)], "1") + + + rc = ReactionSystemWithAutomaton(r,c) + rc.show() + smt_rsc = SmtCheckerRSCParam(rc) + + # Universal property which seems to be true: (holds also existentially) + f1 = Formula_rsLTL.f_G(BagDescription.f_entity("x") > 0, + Formula_rsLTL.f_Implies( + (BagDescription.f_entity('y') == 2), + Formula_rsLTL.f_X( + (BagDescription.f_entity("x") > 1), + BagDescription.f_entity("y") >= 3 + ) + ) + ) + + # lets see if we can find a counterexample to this property: + neg_f1 = Formula_rsLTL.f_F(BagDescription.f_entity("x") > 0, + Formula_rsLTL.f_And( + (BagDescription.f_entity('y') == 2), + Formula_rsLTL.f_X( + (BagDescription.f_entity("x") > 1), + BagDescription.f_entity("y") < 3 + ) + ) + ) + + # we fix the property f1 + # this one holds: + f2 = Formula_rsLTL.f_G(BagDescription.f_entity("x") > 0, + Formula_rsLTL.f_Implies( + (BagDescription.f_entity('y') == 2), + Formula_rsLTL.f_X( + (BagDescription.f_entity("x") > 1) & (BagDescription.f_entity("h") < 1), + BagDescription.f_entity("y") >= 3 + ) + ) + ) + + # neg_f1 = Formula_rsLTL.f_X(BagDescription.f_TRUE(), Formula_rsLTL.f_F(BagDescription.f_entity("x") > 0, + # Formula_rsLTL.f_And( + # (BagDescription.f_entity('y') > 0), + # Formula_rsLTL.f_X( + # BagDescription.f_entity("x") > 0, + # Formula_rsLTL.f_X( + # BagDescription.f_entity("x") > 1, + # BagDescription.f_entity("y") < 3 + # ) + # ) + # ) + # )) + # smt_rsc.check_rsltl(formula=neg_f1) def example44(): 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 index 4604f87..399c3a6 100644 --- a/smt/smt_checker_rsc_param.py +++ b/smt/smt_checker_rsc_param.py @@ -14,14 +14,14 @@ from logics import rsLTL_Encoder # def simplify(x): # return x -class SmtCheckerRSC(object): +class SmtCheckerRSCParam(object): def __init__(self, rsca): - + rsca.sanity_check() - if not rsca.is_with_concentrations(): - raise RuntimeError("RS and CA with concentrations expected") + 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 From c83779c9d5bf26960b502efeee8360abc8318fc9 Mon Sep 17 00:00:00 2001 From: Artur Meski Date: Sun, 13 Aug 2017 15:51:20 +0100 Subject: [PATCH 04/35] Declaration of the intermediate product variables for the reactions and entities that are actually used as products --- ...action_system_with_concentrations_param.py | 22 ++++++--- rs_testing.py | 2 +- smt/smt_checker_rsc_param.py | 48 +++++++++++++++++-- 3 files changed, 62 insertions(+), 10 deletions(-) diff --git a/rs/reaction_system_with_concentrations_param.py b/rs/reaction_system_with_concentrations_param.py index 20e2134..3dde9b3 100644 --- a/rs/reaction_system_with_concentrations_param.py +++ b/rs/reaction_system_with_concentrations_param.py @@ -201,19 +201,29 @@ class ReactionSystemWithConcentrationsParam(ReactionSystem): 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 + # assert False 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)) + producible_entities = self.get_producible_entities() reactions_by_prod = {} diff --git a/rs_testing.py b/rs_testing.py index 4e61f7b..641bf88 100644 --- a/rs_testing.py +++ b/rs_testing.py @@ -90,7 +90,7 @@ def example44_param(): # ) # ) # )) - # smt_rsc.check_rsltl(formula=neg_f1) + smt_rsc.check_rsltl(formula=neg_f1) def example44(): diff --git a/smt/smt_checker_rsc_param.py b/smt/smt_checker_rsc_param.py index 399c3a6..8051744 100644 --- a/smt/smt_checker_rsc_param.py +++ b/smt/smt_checker_rsc_param.py @@ -34,7 +34,13 @@ class SmtCheckerRSCParam(object): self.v = [] self.v_ctx = [] self.ca_state = [] + + # intermediate products: + self.v_improd = [] + self.next_level_to_encode = 0 + + self.producible_entities = self.rs.get_producible_entities() self.loop_position = Int("loop_position") @@ -48,14 +54,15 @@ class SmtCheckerRSCParam(object): self.initialise() def prepare_all_variables(self): - """Encodes all the variables""" + """Prepares all the variables""" self.prepare_state_variables() self.prepare_context_variables() + self.prepare_intermediate_product_variables() self.next_level_to_encode += 1 def prepare_context_variables(self): - """Encodes all the context variables""" + """Prepares all the context variables""" level = self.next_level_to_encode @@ -66,7 +73,7 @@ class SmtCheckerRSCParam(object): self.v_ctx.append(variables) def prepare_state_variables(self): - """Encodes all the state variables""" + """Prepares all the state variables""" level = self.next_level_to_encode @@ -77,6 +84,41 @@ class SmtCheckerRSCParam(object): self.ca_state.append(Int("CA"+str(level)+"_state")) + def prepare_intermediate_product_variables(self): + """ + Prepares the intermediate product variables + carrying the individual concentration levels produced + the reactions. + + These variables are used later on to encode the final + concentration levels for all the entities + """ + + level = self.next_level_to_encode + + if level < 1: + self.v_improd.append([]) + else: + + variables = [] + number_of_reactions = len(self.rs.reactions) + + for reaction in self.rs.reactions: + *_, products = reaction + + reaction_id = self.rs.reactions.index(reaction) + + entities_dict = dict() + for entity, _ in products: + varname = Int("IP" + str(level) + "_R" + + str(reaction_id) + "_e" + str(entity)) + entities_dict[entity] = varname + + variables.append(entities_dict) + + self.v_improd.append(variables) + + def enc_concentration_levels_assertion(self, level): """Encodes assertions that (some) variables need to be >0 From dd4ea5d9396990785c6a8b71a692a47c5e9cfe65 Mon Sep 17 00:00:00 2001 From: Artur Meski Date: Sun, 13 Aug 2017 16:02:22 +0100 Subject: [PATCH 05/35] cleanup --- smt/smt_checker_rsc_param.py | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/smt/smt_checker_rsc_param.py b/smt/smt_checker_rsc_param.py index 8051744..a11e149 100644 --- a/smt/smt_checker_rsc_param.py +++ b/smt/smt_checker_rsc_param.py @@ -28,6 +28,7 @@ class SmtCheckerRSCParam(object): self.initialise() + def initialise(self): """Initialises all the variables used by the checker""" @@ -48,11 +49,13 @@ class SmtCheckerRSCParam(object): self.verification_time = None + def reset(self): """Reinitialises the state of the checker""" self.initialise() + def prepare_all_variables(self): """Prepares all the variables""" @@ -61,6 +64,7 @@ class SmtCheckerRSCParam(object): self.prepare_intermediate_product_variables() self.next_level_to_encode += 1 + def prepare_context_variables(self): """Prepares all the context variables""" @@ -72,6 +76,7 @@ class SmtCheckerRSCParam(object): self.v_ctx.append(variables) + def prepare_state_variables(self): """Prepares all the state variables""" @@ -84,6 +89,7 @@ class SmtCheckerRSCParam(object): self.ca_state.append(Int("CA"+str(level)+"_state")) + def prepare_intermediate_product_variables(self): """ Prepares the intermediate product variables @@ -135,6 +141,7 @@ class SmtCheckerRSCParam(object): 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""" @@ -149,6 +156,7 @@ class SmtCheckerRSCParam(object): 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""" @@ -290,6 +298,7 @@ class SmtCheckerRSCParam(object): 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""" @@ -309,6 +318,7 @@ class SmtCheckerRSCParam(object): return enc_trans + def enc_automaton_trans(self, level): """Encodes the transition relation for the context automaton""" @@ -336,6 +346,7 @@ class SmtCheckerRSCParam(object): return enc_trans + def enc_exact_state(self, level, state): """Encodes the state at the given level with the exact concentration values""" @@ -355,6 +366,7 @@ class SmtCheckerRSCParam(object): # enc = And(enc, self.v[level][entity] == 0) # # return simplify(enc) + def enc_min_state(self, level, state): """Encodes the state at the given level with the minimal required concentration levels""" @@ -370,6 +382,7 @@ class SmtCheckerRSCParam(object): # 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""" @@ -386,6 +399,7 @@ class SmtCheckerRSCParam(object): enc = And(enc, self.v[level][e_id] < conc) return simplify(enc) + def decode_witness(self, max_level, print_model=False): @@ -419,6 +433,7 @@ class SmtCheckerRSCParam(object): 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""" @@ -495,6 +510,7 @@ class SmtCheckerRSCParam(object): 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""" @@ -504,6 +520,7 @@ class SmtCheckerRSCParam(object): 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""" @@ -518,6 +535,7 @@ class SmtCheckerRSCParam(object): eq_enc = simplify(And(eq_enc, eq_enc_ctxaut)) return eq_enc + def get_loop_encodings(self): @@ -537,6 +555,7 @@ class SmtCheckerRSCParam(object): return loop_enc + def check_reachability(self, state, print_witness=True, print_time=True, print_mem=True, max_level=1000): """Main testing function""" @@ -598,9 +617,11 @@ class SmtCheckerRSCParam(object): 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""" @@ -656,3 +677,5 @@ class SmtCheckerRSCParam(object): if not (x == "y" or x == "yes"): break + +# EOF \ No newline at end of file From 9ca4d3954f67b6fed0e575800dc00dc8d1ca9b97 Mon Sep 17 00:00:00 2001 From: Artur Meski Date: Sun, 13 Aug 2017 20:11:25 +0100 Subject: [PATCH 06/35] Initial encoding of reactions --- smt/smt_checker_rsc_param.py | 211 ++++++++++------------------------- 1 file changed, 58 insertions(+), 153 deletions(-) diff --git a/smt/smt_checker_rsc_param.py b/smt/smt_checker_rsc_param.py index a11e149..7c344f9 100644 --- a/smt/smt_checker_rsc_param.py +++ b/smt/smt_checker_rsc_param.py @@ -103,10 +103,15 @@ class SmtCheckerRSCParam(object): level = self.next_level_to_encode if level < 1: - self.v_improd.append([]) + # + # 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.v_improd.append(None) else: - variables = [] + reactions_dict = dict() number_of_reactions = len(self.rs.reactions) for reaction in self.rs.reactions: @@ -120,13 +125,15 @@ class SmtCheckerRSCParam(object): str(reaction_id) + "_e" + str(entity)) entities_dict[entity] = varname - variables.append(entities_dict) - - self.v_improd.append(variables) + reactions_dict[reaction_id] = entities_dict + + self.v_improd.append(reactions_dict) + print(self.v_improd) def enc_concentration_levels_assertion(self, level): - """Encodes assertions that (some) variables need to be >0 + """ + 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. @@ -157,144 +164,7 @@ class SmtCheckerRSCParam(object): 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 == []: - return simplify(self.v[level+1][prod_entity] == 0) # this should never happen - - 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)) - - # for reactants,inhibitors,products in rcts_for_prod_entity: - # enc_reactants = True - # enc_inhibitors = True - # # enc_products -- below - # - # for reactant,concentration in reactants: - # enc_reactants = simplify(And(enc_reactants, - # Or(self.v[level][reactant] >= concentration, self.v_ctx[level][reactant] >= concentration))) - # for inhibitor,concentration in inhibitors: - # enc_inhibitors = simplify(And(enc_inhibitors, - # And(self.v[level][inhibitor] < concentration, self.v_ctx[level][inhibitor] < concentration))) - # - # enc_products = self.v[level+1][products[0][0]] == products[0][1] - # - # enc_enabledness = simplify(Or(enc_enabledness, And(enc_reactants, enc_inhibitors))) - # enc_rct_prod = simplify(Or(enc_rct_prod, And(enc_reactants, enc_inhibitors, enc_products))) - - # -------- 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_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): return simplify(And(self.enc_rs_trans(level), self.enc_automaton_trans(level))) @@ -302,19 +172,54 @@ class SmtCheckerRSCParam(object): def enc_rs_trans(self, level): """Encodes the transition relation""" - unused_entities = set(range(len(self.rs.background_set))) + # + # 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. + # enc_trans = True - reactions = self.rs.get_reactions_by_product() - meta_reactions = self.rs.meta_reactions + for reaction in self.rs.reactions: + + reactants, inhibitors, products = reaction + reaction_id = self.rs.reactions.index(reaction) + + enc_reactants = True + for entity, conc in reactants: + enc_reactants = And(enc_reactants, + Or(self.v[level][entity] >= conc, self.v_ctx[level][entity] >= conc)) + + enc_inhibitors = True + for entity, conc in inhibitors: + enc_inhibitors = And(enc_inhibitors, + And(self.v[level][entity] < conc, self.v_ctx[level][entity] < conc)) + + enc_products = True + # + # + # + for entity, conc in products: + enc_products = And(enc_products, + self.v_improd[level+1][reaction_id][entity] == conc) + + # + # (R and I) iff P + # + enc_reaction = simplify(And(enc_reactants, enc_inhibitors) == enc_products) + + enc_trans = simplify(And(enc_trans, enc_reaction)) + + print(enc_trans) - 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)) + # + # TODO: + # + # Max of all the produced concentrations for each entity/product... + return enc_trans From 15cbf62e8224bfbefc0fa0008c81fae6a255b294 Mon Sep 17 00:00:00 2001 From: Artur Meski Date: Mon, 14 Aug 2017 21:31:43 +0100 Subject: [PATCH 07/35] Max encoding, WIP --- smt/smt_checker_rsc_param.py | 44 ++++++++++++++++++++++++++++++++++-- 1 file changed, 42 insertions(+), 2 deletions(-) diff --git a/smt/smt_checker_rsc_param.py b/smt/smt_checker_rsc_param.py index 7c344f9..7bcd7af 100644 --- a/smt/smt_checker_rsc_param.py +++ b/smt/smt_checker_rsc_param.py @@ -39,6 +39,8 @@ class SmtCheckerRSCParam(object): # intermediate products: self.v_improd = [] + self.v_improd_for_entities = [] + self.next_level_to_encode = 0 self.producible_entities = self.rs.get_producible_entities() @@ -109,27 +111,36 @@ class SmtCheckerRSCParam(object): # which are always at level+1. # self.v_improd.append(None) + self.v_improd_for_entities.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() - for entity, _ in products: + for entity, conc in products: varname = Int("IP" + str(level) + "_R" + str(reaction_id) + "_e" + str(entity)) entities_dict[entity] = varname + all_entities_dict.setdefault(entity, []) + all_entities_dict[entity].append((conc,varname)) + reactions_dict[reaction_id] = entities_dict self.v_improd.append(reactions_dict) + self.v_improd_for_entities.append(all_entities_dict) print(self.v_improd) + + print(self.v_improd_for_entities) def enc_concentration_levels_assertion(self, level): """ @@ -166,7 +177,7 @@ class SmtCheckerRSCParam(object): def enc_transition_relation(self, level): - return simplify(And(self.enc_rs_trans(level), self.enc_automaton_trans(level))) + return simplify(And(self.enc_rs_trans(level), self.enc_automaton_trans(level))) def enc_rs_trans(self, level): @@ -220,9 +231,38 @@ class SmtCheckerRSCParam(object): # # Max of all the produced concentrations for each entity/product... + enc_max_prod = True + + current_v_improd_for_entities = self.v_improd_for_entities[level+1] + for entity, per_reaction_vars in current_v_improd_for_entities.items(): + + # enc_max_single_ent = True + + sorted_vars_by_conc = sorted(per_reaction_vars, key=lambda conc_var: conc_var[0]) + list_of_vars = [v for c,v in sorted_vars_by_conc] + + print(list_of_vars) return enc_trans + # def enc_max(self, elements): + # + # if len(elements) == 1: + # return elements[0] + # + # elif len(elements) > 1: + # + # # If(a > b, a, b) + # def MAX(a, b): + # return If(a > b, a, b) + # + # for i in range(1,len(elements)): + # MAX(elements[]) + # + # else: + # return None + # + def enc_automaton_trans(self, level): """Encodes the transition relation for the context automaton""" From 2b7d5e77dd3a0e226ceddbcdea231f0238ff58d5 Mon Sep 17 00:00:00 2001 From: Artur Meski Date: Tue, 15 Aug 2017 21:41:41 +0100 Subject: [PATCH 08/35] enc_max --- smt/smt_checker_rsc_param.py | 44 +++++++++++++++++++----------------- 1 file changed, 23 insertions(+), 21 deletions(-) diff --git a/smt/smt_checker_rsc_param.py b/smt/smt_checker_rsc_param.py index 7bcd7af..73d412a 100644 --- a/smt/smt_checker_rsc_param.py +++ b/smt/smt_checker_rsc_param.py @@ -14,6 +14,9 @@ from logics import rsLTL_Encoder # def simplify(x): # return x +def MAX(a, b): + return If(a > b, a, b) + class SmtCheckerRSCParam(object): def __init__(self, rsca): @@ -131,7 +134,7 @@ class SmtCheckerRSCParam(object): entities_dict[entity] = varname all_entities_dict.setdefault(entity, []) - all_entities_dict[entity].append((conc,varname)) + all_entities_dict[entity].append(varname) reactions_dict[reaction_id] = entities_dict @@ -238,31 +241,30 @@ class SmtCheckerRSCParam(object): # enc_max_single_ent = True - sorted_vars_by_conc = sorted(per_reaction_vars, key=lambda conc_var: conc_var[0]) - list_of_vars = [v for c,v in sorted_vars_by_conc] + #sorted_vars_by_conc = sorted(per_reaction_vars, key=lambda conc_var: conc_var[0]) + #list_of_vars = [v for c,v in sorted_vars_by_conc] - print(list_of_vars) + print(per_reaction_vars, "--->", self.enc_max(per_reaction_vars)) return enc_trans - # def enc_max(self, elements): - # - # if len(elements) == 1: - # return elements[0] - # - # elif len(elements) > 1: - # - # # If(a > b, a, b) - # def MAX(a, b): - # return If(a > b, a, b) - # - # for i in range(1,len(elements)): - # MAX(elements[]) - # - # else: - # return None - # + def enc_max(self, elements): + + if len(elements) == 1: + return MAX(0, elements[0]) + + elif len(elements) > 1: + + enc = 0 + for i in range(len(elements) - 1): + enc = MAX(enc, MAX(elements[i], elements[i + 1])) + + return enc + + else: + return None + def enc_automaton_trans(self, level): """Encodes the transition relation for the context automaton""" From ce6a1f6f195c0a50965be22334c9eea1ba7a1250 Mon Sep 17 00:00:00 2001 From: Artur Meski Date: Sun, 20 Aug 2017 19:57:23 +0100 Subject: [PATCH 09/35] transition relation with max calculation --- ...action_system_with_concentrations_param.py | 3 +- smt/smt_checker_rsc_param.py | 90 ++++++++++--------- 2 files changed, 50 insertions(+), 43 deletions(-) diff --git a/rs/reaction_system_with_concentrations_param.py b/rs/reaction_system_with_concentrations_param.py index 3dde9b3..108bb12 100644 --- a/rs/reaction_system_with_concentrations_param.py +++ b/rs/reaction_system_with_concentrations_param.py @@ -64,7 +64,8 @@ class ReactionSystemWithConcentrationsParam(ReactionSystem): def get_state_ids(self, state): """Returns entities of the given state without levels""" - return [e for e,c in state] + # 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: diff --git a/smt/smt_checker_rsc_param.py b/smt/smt_checker_rsc_param.py index 73d412a..a71a0ae 100644 --- a/smt/smt_checker_rsc_param.py +++ b/smt/smt_checker_rsc_param.py @@ -47,6 +47,7 @@ class SmtCheckerRSCParam(object): self.next_level_to_encode = 0 self.producible_entities = self.rs.get_producible_entities() + self.improducible_entities = set(self.rs.get_state_ids(self.rs.background_set)) - self.producible_entities self.loop_position = Int("loop_position") @@ -141,9 +142,9 @@ class SmtCheckerRSCParam(object): self.v_improd.append(reactions_dict) self.v_improd_for_entities.append(all_entities_dict) - print(self.v_improd) + # print(self.v_improd) - print(self.v_improd_for_entities) + # print(self.v_improd_for_entities) def enc_concentration_levels_assertion(self, level): """ @@ -162,26 +163,26 @@ class SmtCheckerRSCParam(object): 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]: - rs_init_state_enc = simplify(And(rs_init_state_enc, v == 0)) # the initial concentration levels are zeroed + # 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_transition_relation(self, level): - return simplify(And(self.enc_rs_trans(level), self.enc_automaton_trans(level))) + 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""" @@ -191,68 +192,76 @@ class SmtCheckerRSCParam(object): # # 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. # enc_trans = True for reaction in self.rs.reactions: - + reactants, inhibitors, products = reaction reaction_id = self.rs.reactions.index(reaction) - + enc_reactants = True for entity, conc in reactants: - enc_reactants = And(enc_reactants, - Or(self.v[level][entity] >= conc, self.v_ctx[level][entity] >= conc)) - + enc_reactants = And(enc_reactants, Or( + self.v[level][entity] >= conc, self.v_ctx[level][entity] >= conc)) + enc_inhibitors = True for entity, conc in inhibitors: - enc_inhibitors = And(enc_inhibitors, - And(self.v[level][entity] < conc, self.v_ctx[level][entity] < conc)) - + enc_inhibitors = And(enc_inhibitors, And( + self.v[level][entity] < conc, self.v_ctx[level][entity] < conc)) + enc_products = True - # - # - # for entity, conc in products: - enc_products = And(enc_products, - self.v_improd[level+1][reaction_id][entity] == conc) - + enc_products = And(enc_products, self.v_improd[ + level + 1][reaction_id][entity] == conc) + # # (R and I) iff P # - enc_reaction = simplify(And(enc_reactants, enc_inhibitors) == enc_products) - + enc_reaction = simplify( + And(enc_reactants, enc_inhibitors) == enc_products) + enc_trans = simplify(And(enc_trans, enc_reaction)) - - print(enc_trans) + + # print(enc_trans) # # TODO: # # Max of all the produced concentrations for each entity/product... - + enc_max_prod = True - current_v_improd_for_entities = self.v_improd_for_entities[level+1] + current_v_improd_for_entities = self.v_improd_for_entities[level + 1] for entity, per_reaction_vars in current_v_improd_for_entities.items(): # enc_max_single_ent = True - #sorted_vars_by_conc = sorted(per_reaction_vars, key=lambda conc_var: conc_var[0]) - #list_of_vars = [v for c,v in sorted_vars_by_conc] + # sorted_vars_by_conc = sorted(per_reaction_vars, key=lambda conc_var: conc_var[0]) + # list_of_vars = [v for c,v in sorted_vars_by_conc] - print(per_reaction_vars, "--->", self.enc_max(per_reaction_vars)) + # print(per_reaction_vars, "--->", self.enc_max(per_reaction_vars)) + + enc_max_prod = simplify( + And(enc_max_prod, self.v[level + 1][entity] == self.enc_max(per_reaction_vars))) + + for entity in self.improducible_entities: + enc_max_prod = simplify( + And(enc_max_prod, self.v[level + 1][entity] == 0)) + + enc_trans_with_max = simplify(And(enc_max_prod, enc_trans)) + + return enc_trans_with_max - return enc_trans - - def enc_max(self, elements): + enc = None + if len(elements) == 1: - return MAX(0, elements[0]) + enc = MAX(0, elements[0]) elif len(elements) > 1: @@ -260,11 +269,8 @@ class SmtCheckerRSCParam(object): for i in range(len(elements) - 1): enc = MAX(enc, MAX(elements[i], elements[i + 1])) - return enc + return enc - else: - return None - def enc_automaton_trans(self, level): """Encodes the transition relation for the context automaton""" From 203182c520b732ae9a984d820ef41af573aae07d Mon Sep 17 00:00:00 2001 From: Artur Meski Date: Sun, 3 Sep 2017 14:53:08 +0100 Subject: [PATCH 10/35] parametric example --- rs_testing.py | 272 ++++++++++++++++++++++++++++---------------------- 1 file changed, 150 insertions(+), 122 deletions(-) diff --git a/rs_testing.py b/rs_testing.py index 641bf88..6bb523d 100644 --- a/rs_testing.py +++ b/rs_testing.py @@ -12,150 +12,101 @@ def run_tests(): # heat_shock_response() # scalable_chain(print_system=True) # example44() - example44_param() + # example44_param() + trivial_param() + + +def trivial_param(): -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)]) - + 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)) + + r.add_reaction([("x", 1)], [("c", 1)], [("y", 2)]) + r.add_reaction([("y", 1)], [("c", 1)], [("z", 1)]) + r.add_reaction([("z", 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 = ReactionSystemWithAutomaton(r, c) rc.show() smt_rsc = SmtCheckerRSCParam(rc) - - # Universal property which seems to be true: (holds also existentially) - f1 = Formula_rsLTL.f_G(BagDescription.f_entity("x") > 0, - Formula_rsLTL.f_Implies( - (BagDescription.f_entity('y') == 2), - Formula_rsLTL.f_X( - (BagDescription.f_entity("x") > 1), - BagDescription.f_entity("y") >= 3 - ) - ) - ) - - # lets see if we can find a counterexample to this property: - neg_f1 = Formula_rsLTL.f_F(BagDescription.f_entity("x") > 0, - Formula_rsLTL.f_And( - (BagDescription.f_entity('y') == 2), - Formula_rsLTL.f_X( - (BagDescription.f_entity("x") > 1), - BagDescription.f_entity("y") < 3 - ) - ) - ) - - # we fix the property f1 - # this one holds: - f2 = Formula_rsLTL.f_G(BagDescription.f_entity("x") > 0, - Formula_rsLTL.f_Implies( - (BagDescription.f_entity('y') == 2), - Formula_rsLTL.f_X( - (BagDescription.f_entity("x") > 1) & (BagDescription.f_entity("h") < 1), - BagDescription.f_entity("y") >= 3 - ) - ) - ) - - # neg_f1 = Formula_rsLTL.f_X(BagDescription.f_TRUE(), Formula_rsLTL.f_F(BagDescription.f_entity("x") > 0, - # Formula_rsLTL.f_And( - # (BagDescription.f_entity('y') > 0), - # Formula_rsLTL.f_X( - # BagDescription.f_entity("x") > 0, - # Formula_rsLTL.f_X( - # BagDescription.f_entity("x") > 1, - # BagDescription.f_entity("y") < 3 - # ) - # ) - # ) - # )) - smt_rsc.check_rsltl(formula=neg_f1) - -def example44(): - - 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)]) - + + f2 = Formula_rsLTL.f_F( + BagDescription.f_TRUE(), + BagDescription.f_entity("final") >= 1) + + smt_rsc.check_rsltl(formula=f2) + + +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") + 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 = 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), @@ -169,6 +120,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): From 4b9aa59affd24cc31157bde051e3699249a993e7 Mon Sep 17 00:00:00 2001 From: Artur Meski Date: Sun, 3 Sep 2017 14:53:39 +0100 Subject: [PATCH 11/35] some sanity checks for a common mistake --- logics/bags.py | 60 ++++++++++++++++++++++------------------- logics/rsltl.py | 71 ++++++++++++++++++++++++++++++------------------- 2 files changed, 77 insertions(+), 54 deletions(-) diff --git a/logics/bags.py b/logics/bags.py index 30ef358..b69d694 100644 --- a/logics/bags.py +++ b/logics/bags.py @@ -1,15 +1,18 @@ 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 @@ -31,45 +34,48 @@ 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: " + str(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) - + 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..83d61d1 100644 --- a/logics/rsltl.py +++ b/logics/rsltl.py @@ -3,22 +3,30 @@ 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 - + + 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) - + def __repr__(self): if self.f_type == rsLTL_form_type.bag: return repr(self.bag_descr) @@ -37,9 +45,15 @@ class Formula_rsLTL(object): if self.f_type == rsLTL_form_type.l_implies: return "( " + repr(self.left_operand) + " => " + repr(self.right_operand) + " )" if self.f_type == rsLTL_form_type.t_until: - return "( " + repr(self.left_operand) + " U[" + repr(self.sub_operand) + "] " + repr(self.right_operand) + " )" + return "( " + repr( + self.left_operand) + " U[" + repr( + self.sub_operand) + "] " + repr( + self.right_operand) + " )" if self.f_type == rsLTL_form_type.t_release: - return "( " + repr(self.left_operand) + " R[" + repr(self.sub_operand) + "] " + repr(self.right_operand) + " )" + return "( " + repr( + self.left_operand) + " R[" + repr( + self.sub_operand) + "] " + repr( + self.right_operand) + " )" @property def is_bag(self): @@ -47,48 +61,51 @@ 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_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 From ac150b2744a07048b29353d41705876fad1dcabc Mon Sep 17 00:00:00 2001 From: Artur Meski Date: Sun, 3 Sep 2017 14:54:10 +0100 Subject: [PATCH 12/35] cleanup --- rs/reaction_system_with_concentrations_param.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/rs/reaction_system_with_concentrations_param.py b/rs/reaction_system_with_concentrations_param.py index 108bb12..405387e 100644 --- a/rs/reaction_system_with_concentrations_param.py +++ b/rs/reaction_system_with_concentrations_param.py @@ -3,6 +3,10 @@ from colour import * from rs.reaction_system import ReactionSystem +class ParameterSet(object): + def __init__(self, name): + self.name = name + class ReactionSystemWithConcentrationsParam(ReactionSystem): def __init__(self): From a2dc2ec31fdcf969040f12d9dc96a9827a1bf1e4 Mon Sep 17 00:00:00 2001 From: Artur Meski Date: Sun, 3 Sep 2017 15:13:23 +0100 Subject: [PATCH 13/35] cleanup --- ...action_system_with_concentrations_param.py | 52 ++++++++++--------- 1 file changed, 28 insertions(+), 24 deletions(-) diff --git a/rs/reaction_system_with_concentrations_param.py b/rs/reaction_system_with_concentrations_param.py index 405387e..bf9d7a8 100644 --- a/rs/reaction_system_with_concentrations_param.py +++ b/rs/reaction_system_with_concentrations_param.py @@ -4,32 +4,35 @@ from colour import * from rs.reaction_system import ReactionSystem class ParameterSet(object): + def __init__(self, name): self.name = name + class ReactionSystemWithConcentrationsParam(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() + 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 + 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") + raise RuntimeError( + "Bad entity type when adding background set element") self.assume_not_in_bgset(name) self.background_set.append(name) @@ -43,7 +46,7 @@ class ReactionSystemWithConcentrationsParam(ReactionSystem): 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: @@ -51,11 +54,11 @@ class ReactionSystemWithConcentrationsParam(ReactionSystem): 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 @@ -63,7 +66,7 @@ class ReactionSystemWithConcentrationsParam(ReactionSystem): print("FATAL. Invalid entity+concentration:") print(e) exit(1) - + return False def get_state_ids(self, state): @@ -73,38 +76,39 @@ class ReactionSystemWithConcentrationsParam(ReactionSystem): def has_non_zero_concentration(self, elem): if elem[1] < 1: - raise RuntimeError("Unexpected concentration level in state: " + str(elem)) + 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)) + 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)) + 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 + 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 From 2957a7ac520c5b0084c824f935d92a3e3a5c73e9 Mon Sep 17 00:00:00 2001 From: Artur Meski Date: Sun, 3 Sep 2017 15:15:26 +0100 Subject: [PATCH 14/35] cleanup --- ...action_system_with_concentrations_param.py | 121 ++++++++++-------- 1 file changed, 67 insertions(+), 54 deletions(-) diff --git a/rs/reaction_system_with_concentrations_param.py b/rs/reaction_system_with_concentrations_param.py index bf9d7a8..328afbb 100644 --- a/rs/reaction_system_with_concentrations_param.py +++ b/rs/reaction_system_with_concentrations_param.py @@ -274,47 +274,47 @@ class ReactionSystemWithConcentrationsParam(ReactionSystem): """ Translates RSC into RS """ - + rs = ReactionSystem() - - for reactants,inhibitors,products in self.reactions: + + for reactants, inhibitors, products in self.reactions: new_reactants = [] new_inhibitors = [] new_products = [] - for ent,conc in reactants: + 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: + + 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): + 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: - + + 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: + + 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: + + for ent, conc in inhibitors: n = self.get_entity_name(ent) + "#" + str(conc) rs.ensure_bg_set_entity(n) new_inhibitors.append(n) @@ -323,81 +323,94 @@ class ReactionSystemWithConcentrationsParam(ReactionSystem): 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( + "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): + + 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): + + 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 + 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)) + 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): + 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 + 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)) - + 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(): - + 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( + "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): - + + 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: + + 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) + new_inhibitors.append(n) - for i in range(1,value+1): + for i in range(1, value + 1): new_products.append(e_value(i)) - - rs.add_reaction(new_reactants,new_inhibitors,new_products) - + + rs.add_reaction(new_reactants, new_inhibitors, new_products) + return rs - + # class ReactionSystemWithAutomaton(object): # From 782ae117a3ac1a88b2bb883685046f55be145375 Mon Sep 17 00:00:00 2001 From: Artur Meski Date: Sun, 3 Sep 2017 17:59:26 +0100 Subject: [PATCH 15/35] introducing parameters --- rs/__init__.py | 4 +- ...action_system_with_concentrations_param.py | 90 +++++++++++++------ rs_testing.py | 4 +- smt/smt_checker_rsc_param.py | 20 ++--- 4 files changed, 75 insertions(+), 43 deletions(-) diff --git a/rs/__init__.py b/rs/__init__.py index 01001bf..2e56fbd 100644 --- a/rs/__init__.py +++ b/rs/__init__.py @@ -2,11 +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 +from rs.reaction_system_with_concentrations_param import ReactionSystemWithConcentrationsParam, ParameterObj 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_concentrations_param.py b/rs/reaction_system_with_concentrations_param.py index 328afbb..c86194f 100644 --- a/rs/reaction_system_with_concentrations_param.py +++ b/rs/reaction_system_with_concentrations_param.py @@ -3,10 +3,13 @@ from colour import * from rs.reaction_system import ReactionSystem -class ParameterSet(object): +class ParameterObj(object): def __init__(self, name): self.name = name + + def __repr__(self): + return "@{0}".format(self.name) class ReactionSystemWithConcentrationsParam(ReactionSystem): @@ -14,6 +17,7 @@ class ReactionSystemWithConcentrationsParam(ReactionSystem): def __init__(self): self.reactions = [] + self.parametric_reactions = [] self.meta_reactions = dict() self.permanent_entities = dict() self.background_set = [] @@ -85,31 +89,55 @@ class ReactionSystemWithConcentrationsParam(ReactionSystem): if R == [] and not ignore_empty_R: raise RuntimeError("No reactants defined") + # + # REACTANTS + # 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 + 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 = [] - 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 + 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 = [] - 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)) + 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 @@ -119,7 +147,11 @@ class ReactionSystemWithConcentrationsParam(ReactionSystem): if P == []: raise RuntimeError("No products defined") reaction = self.process_rip(R,I,P) - self.reactions.append(reaction) + + if self.is_parametric_reaction(reaction): + self.parametric_reactions.append(reaction) + else: + self.reactions.append(reaction) def add_reaction_without_reactants(self, R, I, P): """Adds a reaction""" @@ -174,10 +206,18 @@ class ReactionSystemWithConcentrationsParam(ReactionSystem): 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] + """ + 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): diff --git a/rs_testing.py b/rs_testing.py index 6bb523d..5b77c7c 100644 --- a/rs_testing.py +++ b/rs_testing.py @@ -27,7 +27,7 @@ def trivial_param(): r.add_reaction([("x", 1)], [("c", 1)], [("y", 2)]) r.add_reaction([("y", 1)], [("c", 1)], [("z", 1)]) - r.add_reaction([("z", 1)], [("c", 1)], [("final", 1)]) + r.add_reaction(ParameterObj("P1"), [("c", 1)], [("final", 1)]) c = ContextAutomatonWithConcentrations(r) c.add_init_state("0") @@ -43,7 +43,7 @@ def trivial_param(): BagDescription.f_TRUE(), BagDescription.f_entity("final") >= 1) - smt_rsc.check_rsltl(formula=f2) + smt_rsc.check_rsltl(formula=f2, max_level=5) def example44_param(): diff --git a/smt/smt_checker_rsc_param.py b/smt/smt_checker_rsc_param.py index a71a0ae..50ae1bb 100644 --- a/smt/smt_checker_rsc_param.py +++ b/smt/smt_checker_rsc_param.py @@ -14,7 +14,7 @@ from logics import rsLTL_Encoder # def simplify(x): # return x -def MAX(a, b): +def z3_max(a, b): return If(a > b, a, b) class SmtCheckerRSCParam(object): @@ -31,7 +31,6 @@ class SmtCheckerRSCParam(object): self.initialise() - def initialise(self): """Initialises all the variables used by the checker""" @@ -41,7 +40,6 @@ class SmtCheckerRSCParam(object): # intermediate products: self.v_improd = [] - self.v_improd_for_entities = [] self.next_level_to_encode = 0 @@ -55,12 +53,10 @@ class SmtCheckerRSCParam(object): self.verification_time = None - def reset(self): """Reinitialises the state of the checker""" - self.initialise() - + self.initialise() def prepare_all_variables(self): """Prepares all the variables""" @@ -68,8 +64,7 @@ class SmtCheckerRSCParam(object): self.prepare_state_variables() self.prepare_context_variables() self.prepare_intermediate_product_variables() - self.next_level_to_encode += 1 - + self.next_level_to_encode += 1 def prepare_context_variables(self): """Prepares all the context variables""" @@ -82,7 +77,6 @@ class SmtCheckerRSCParam(object): self.v_ctx.append(variables) - def prepare_state_variables(self): """Prepares all the state variables""" @@ -95,7 +89,6 @@ class SmtCheckerRSCParam(object): self.ca_state.append(Int("CA"+str(level)+"_state")) - def prepare_intermediate_product_variables(self): """ Prepares the intermediate product variables @@ -261,16 +254,15 @@ class SmtCheckerRSCParam(object): enc = None if len(elements) == 1: - enc = MAX(0, elements[0]) + enc = z3_max(0, elements[0]) elif len(elements) > 1: enc = 0 for i in range(len(elements) - 1): - enc = MAX(enc, MAX(elements[i], elements[i + 1])) + enc = z3_max(enc, z3_max(elements[i], elements[i + 1])) return enc - def enc_automaton_trans(self, level): """Encodes the transition relation for the context automaton""" @@ -603,7 +595,7 @@ class SmtCheckerRSCParam(object): 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))) From 3cba5901f9d0aa608823c0197738e5a07b9fbde9 Mon Sep 17 00:00:00 2001 From: Artur Meski Date: Sun, 3 Sep 2017 18:52:25 +0100 Subject: [PATCH 16/35] get_param --- ...action_system_with_concentrations_param.py | 19 +++++++++++++++++++ rs_testing.py | 4 +++- 2 files changed, 22 insertions(+), 1 deletion(-) diff --git a/rs/reaction_system_with_concentrations_param.py b/rs/reaction_system_with_concentrations_param.py index c86194f..f984bab 100644 --- a/rs/reaction_system_with_concentrations_param.py +++ b/rs/reaction_system_with_concentrations_param.py @@ -17,6 +17,7 @@ class ReactionSystemWithConcentrationsParam(ReactionSystem): def __init__(self): self.reactions = [] + self.parameters = dict() self.parametric_reactions = [] self.meta_reactions = dict() self.permanent_entities = dict() @@ -49,6 +50,24 @@ class ReactionSystemWithConcentrationsParam(ReactionSystem): 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: diff --git a/rs_testing.py b/rs_testing.py index 5b77c7c..effc5a3 100644 --- a/rs_testing.py +++ b/rs_testing.py @@ -24,10 +24,12 @@ def trivial_param(): 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(ParameterObj("P1"), [("c", 1)], [("final", 1)]) + r.add_reaction(param_p1, [("c", 1)], [("final", 1)]) c = ContextAutomatonWithConcentrations(r) c.add_init_state("0") From 690d3e5d0a83635c86ecf3b18458998608d24781 Mon Sep 17 00:00:00 2001 From: Artur Meski Date: Sun, 3 Sep 2017 19:07:34 +0100 Subject: [PATCH 17/35] comment --- rs_testing.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/rs_testing.py b/rs_testing.py index effc5a3..9f9110c 100644 --- a/rs_testing.py +++ b/rs_testing.py @@ -41,11 +41,14 @@ def trivial_param(): rc.show() smt_rsc = SmtCheckerRSCParam(rc) - f2 = Formula_rsLTL.f_F( + f1 = Formula_rsLTL.f_F( BagDescription.f_TRUE(), BagDescription.f_entity("final") >= 1) - - smt_rsc.check_rsltl(formula=f2, max_level=5) + + # + # WARNING: depth limit is set + # + smt_rsc.check_rsltl(formula=f1, max_level=5) def example44_param(): From 72dfd3b41a511adb1f584eca5831a517fa5c6e71 Mon Sep 17 00:00:00 2001 From: Artur Meski Date: Sun, 10 Sep 2017 18:22:42 +0100 Subject: [PATCH 18/35] fixed permissions --- rs_examples.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) mode change 100755 => 100644 rs_examples.py diff --git a/rs_examples.py b/rs_examples.py old mode 100755 new mode 100644 From dab97e9daefeec7721b37757a78c5eafcfac0cbf Mon Sep 17 00:00:00 2001 From: Artur Meski Date: Sun, 10 Sep 2017 21:47:57 +0100 Subject: [PATCH 19/35] parameter variables; some cleanup --- colour.py | 4 +- rs_examples.py | 1 + rs_testing.py | 2 +- smt/smt_checker_rsc_param.py | 144 +++++++++++++++++++++++------------ 4 files changed, 102 insertions(+), 49 deletions(-) mode change 100644 => 100755 rs_examples.py diff --git a/colour.py b/colour.py index 07cc048..29ba737 100644 --- a/colour.py +++ b/colour.py @@ -14,4 +14,6 @@ def colour_str(col, s): return col + 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) + +# EOF diff --git a/rs_examples.py b/rs_examples.py old mode 100644 new mode 100755 index bbe9478..f0e26b9 --- 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 9f9110c..4d746c6 100644 --- a/rs_testing.py +++ b/rs_testing.py @@ -48,7 +48,7 @@ def trivial_param(): # # WARNING: depth limit is set # - smt_rsc.check_rsltl(formula=f1, max_level=5) + smt_rsc.check_rsltl(formula=f1, max_level=10, print_witness=True) def example44_param(): diff --git a/smt/smt_checker_rsc_param.py b/smt/smt_checker_rsc_param.py index 50ae1bb..fa9941f 100644 --- a/smt/smt_checker_rsc_param.py +++ b/smt/smt_checker_rsc_param.py @@ -41,17 +41,25 @@ class SmtCheckerRSCParam(object): # intermediate products: self.v_improd = [] self.v_improd_for_entities = [] + + self.v_param = [] 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 + # 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) self.loop_position = Int("loop_position") self.solver = Solver() #For("QF_FD") - self.verification_time = None + self.verification_time = None def reset(self): """Reinitialises the state of the checker""" @@ -64,6 +72,7 @@ class SmtCheckerRSCParam(object): self.prepare_state_variables() self.prepare_context_variables() self.prepare_intermediate_product_variables() + self.prepare_param_variables() self.next_level_to_encode += 1 def prepare_context_variables(self): @@ -93,7 +102,7 @@ class SmtCheckerRSCParam(object): """ Prepares the intermediate product variables carrying the individual concentration levels produced - the reactions. + by the reactions. These variables are used later on to encode the final concentration levels for all the entities @@ -123,12 +132,12 @@ class SmtCheckerRSCParam(object): entities_dict = dict() for entity, conc in products: - varname = Int("IP" + str(level) + "_R" + + new_var = Int("IP" + str(level) + "_R" + str(reaction_id) + "_e" + str(entity)) - entities_dict[entity] = varname + entities_dict[entity] = new_var all_entities_dict.setdefault(entity, []) - all_entities_dict[entity].append(varname) + all_entities_dict[entity].append(new_var) reactions_dict[reaction_id] = entities_dict @@ -139,6 +148,25 @@ class SmtCheckerRSCParam(object): # print(self.v_improd_for_entities) + 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. + """ + + level = self.next_level_to_encode + + params = [] + for entity in self.rs.background_set: + new_var = Int("L{:d}_Pm_{:s}".format(level, entity)) + params.append(new_var) + self.v_param.append(params) + + print(self.v_param) + def enc_concentration_levels_assertion(self, level): """ Encodes assertions that (some) variables need to be >0 @@ -177,6 +205,46 @@ class SmtCheckerRSCParam(object): And(self.enc_rs_trans(level), self.enc_automaton_trans(level))) + def enc_single_reaction(self, level, 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) + + enc_reactants = True + for entity, conc in reactants: + enc_reactants = And(enc_reactants, Or( + self.v[level][entity] >= conc, self.v_ctx[level][entity] >= conc)) + + enc_inhibitors = True + for entity, conc in inhibitors: + enc_inhibitors = And(enc_inhibitors, And( + self.v[level][entity] < conc, self.v_ctx[level][entity] < conc)) + + enc_products = True + for entity, conc in products: + enc_products = And(enc_products, self.v_improd[ + level + 1][reaction_id][entity] == conc) + + # + # (R and I) iff P + # + enc_reaction = simplify( + And(enc_reactants, enc_inhibitors) == enc_products) + + return enc_reaction + def enc_rs_trans(self, level): """Encodes the transition relation""" @@ -188,62 +256,43 @@ class SmtCheckerRSCParam(object): # # 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: - - reactants, inhibitors, products = reaction - reaction_id = self.rs.reactions.index(reaction) - - enc_reactants = True - for entity, conc in reactants: - enc_reactants = And(enc_reactants, Or( - self.v[level][entity] >= conc, self.v_ctx[level][entity] >= conc)) - - enc_inhibitors = True - for entity, conc in inhibitors: - enc_inhibitors = And(enc_inhibitors, And( - self.v[level][entity] < conc, self.v_ctx[level][entity] < conc)) - - enc_products = True - for entity, conc in products: - enc_products = And(enc_products, self.v_improd[ - level + 1][reaction_id][entity] == conc) - - # - # (R and I) iff P - # - enc_reaction = simplify( - And(enc_reactants, enc_inhibitors) == enc_products) - + enc_reaction = self.enc_single_reaction(level, reaction) enc_trans = simplify(And(enc_trans, enc_reaction)) - # print(enc_trans) - - # - # TODO: - # - # Max of all the produced concentrations for each entity/product... + # 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.v_improd_for_entities[level + 1] + for entity, per_reaction_vars in current_v_improd_for_entities.items(): - # enc_max_single_ent = True - - # sorted_vars_by_conc = sorted(per_reaction_vars, key=lambda conc_var: conc_var[0]) - # list_of_vars = [v for c,v in sorted_vars_by_conc] - - # print(per_reaction_vars, "--->", self.enc_max(per_reaction_vars)) - enc_max_prod = simplify( And(enc_max_prod, self.v[level + 1][entity] == self.enc_max(per_reaction_vars))) - for entity in self.improducible_entities: - enc_max_prod = simplify( - And(enc_max_prod, self.v[level + 1][entity] == 0)) + # for entity in self.improducible_entities: + # enc_max_prod = simplify( + # And(enc_max_prod, self.v[level + 1][entity] == 0)) enc_trans_with_max = simplify(And(enc_max_prod, enc_trans)) @@ -428,6 +477,7 @@ class SmtCheckerRSCParam(object): result = self.solver.check() if result == sat: print("[" + colour_str(C_BOLD, "+") + "] " + colour_str(C_GREEN, "SAT at level=" + str(self.current_level))) + print(self.solver.model()) if print_witness: print("\n{:=^70}".format("[ WITNESS ]")) self.decode_witness(self.current_level) From 9abd631e4fae806c097cb941711780702739739b Mon Sep 17 00:00:00 2001 From: Artur Meski Date: Sun, 10 Sep 2017 22:02:31 +0100 Subject: [PATCH 20/35] formatting --- smt/smt_checker_rsc_param.py | 97 +++++++++++++++++++++++------------- 1 file changed, 63 insertions(+), 34 deletions(-) diff --git a/smt/smt_checker_rsc_param.py b/smt/smt_checker_rsc_param.py index fa9941f..e00f5b8 100644 --- a/smt/smt_checker_rsc_param.py +++ b/smt/smt_checker_rsc_param.py @@ -402,7 +402,7 @@ class SmtCheckerRSCParam(object): if print_model: print(m) - for level in range(max_level+1): + for level in range(max_level + 1): print("\n{: >70}".format("[ level=" + repr(level) + " ]")) @@ -410,9 +410,12 @@ class SmtCheckerRSCParam(object): 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") + 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( + " " + self.rs.get_entity_name(var_id) + "=" + var_rep, + end="") # print(" " + repr(m[self.v[level][var_id]]), end="") print(" }") @@ -422,20 +425,26 @@ class SmtCheckerRSCParam(object): 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") + 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( + " " + 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): + 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") + + "] 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 @@ -445,38 +454,49 @@ class SmtCheckerRSCParam(object): self.current_level = 0 self.prepare_all_variables() - + self.solver.add(self.enc_concentration_levels_assertion(0)) - + encoder = rsLTL_Encoder(self) - + 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) + " ]")) + 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...") - + + 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...") - + + 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...") + + 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))) + print( + "[" + colour_str(C_BOLD, "+") + "] " + + colour_str( + C_GREEN, "SAT at level=" + str(self.current_level))) print(self.solver.model()) if print_witness: print("\n{:=^70}".format("[ WITNESS ]")) @@ -485,10 +505,12 @@ class SmtCheckerRSCParam(object): else: self.solver.pop() - print("[" + colour_str(C_BOLD, "i") + "] Unrolling the transition relation") + 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 ]")) + 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: @@ -498,13 +520,20 @@ class SmtCheckerRSCParam(object): if print_time: # stop = time() stop = resource.getrusage(resource.RUSAGE_SELF).ru_utime - self.verification_time = stop-start + self.verification_time = stop - start print() - print("\n[i] {: >60}".format(" Time: " + repr(self.verification_time) + " s")) - + 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")) - + 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""" From 58d65453b6c4751cc6d397aca74f089f48aa228e Mon Sep 17 00:00:00 2001 From: Artur Meski Date: Sun, 10 Sep 2017 22:46:06 +0100 Subject: [PATCH 21/35] comments, etc. --- smt/smt_checker_rsc_param.py | 33 +++++++++++++++++++++++++++++---- 1 file changed, 29 insertions(+), 4 deletions(-) diff --git a/smt/smt_checker_rsc_param.py b/smt/smt_checker_rsc_param.py index e00f5b8..b0ec3bb 100644 --- a/smt/smt_checker_rsc_param.py +++ b/smt/smt_checker_rsc_param.py @@ -46,8 +46,8 @@ class SmtCheckerRSCParam(object): self.next_level_to_encode = 0 - # this is probably not needed anymore - self.producible_entities = self.rs.get_producible_entities() + # 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 @@ -156,7 +156,13 @@ class SmtCheckerRSCParam(object): therefore we need separate variables for each element of the background set. """ - + + # TODO: + # + # We should be checking here if we acutally need + # any variables for parameters, because there might + # be no parameters in RS + level = self.next_level_to_encode params = [] @@ -208,7 +214,7 @@ class SmtCheckerRSCParam(object): def enc_single_reaction(self, level, reaction): """ Encodes a single reaction - + For encoding the products we use intermediate variables: * each reaction has its own product variables, @@ -222,20 +228,39 @@ class SmtCheckerRSCParam(object): # we need reaction_id to find the intermediate product variable reaction_id = self.rs.reactions.index(reaction) + # ** REACTANTS ** + enc_reactants = True for entity, conc in reactants: enc_reactants = And(enc_reactants, Or( self.v[level][entity] >= conc, self.v_ctx[level][entity] >= conc)) + # ** INHIBITORS ** + enc_inhibitors = True for entity, conc in inhibitors: enc_inhibitors = And(enc_inhibitors, And( self.v[level][entity] < conc, self.v_ctx[level][entity] < conc)) + # ** PRODUCTS ** + + # entities that were produced + produced = [] enc_products = True for entity, conc in products: enc_products = And(enc_products, self.v_improd[ level + 1][reaction_id][entity] == conc) + produced.append(entity) + + print("----", str(produced)) + + # encoding of the zeros (for the products) + # we iterate through all the entities' indices + + # for entity in range(len(self.rs.background_set)): + # if entity not in produced: + # enc_products = And(enc_products, self.v_improd[ + # level + 1][reaction_id][entity] == 0) # # (R and I) iff P From f7cfeb43c113440ccb853e3fa1696e8430ef2bb3 Mon Sep 17 00:00:00 2001 From: Artur Meski Date: Thu, 14 Sep 2017 22:25:21 +0100 Subject: [PATCH 22/35] parameters - z3... variables --- rs/__init__.py | 2 +- ...action_system_with_concentrations_param.py | 18 ++++++- smt/smt_checker_rsc_param.py | 49 ++++++++++++++----- 3 files changed, 55 insertions(+), 14 deletions(-) diff --git a/rs/__init__.py b/rs/__init__.py index 2e56fbd..129a278 100644 --- a/rs/__init__.py +++ b/rs/__init__.py @@ -2,7 +2,7 @@ from rs.reaction_system import ReactionSystem from rs.context_automaton import ContextAutomaton from rs.reaction_system_with_concentrations import ReactionSystemWithConcentrations -from rs.reaction_system_with_concentrations_param import ReactionSystemWithConcentrationsParam, ParameterObj +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 diff --git a/rs/reaction_system_with_concentrations_param.py b/rs/reaction_system_with_concentrations_param.py index f984bab..c06bdc3 100644 --- a/rs/reaction_system_with_concentrations_param.py +++ b/rs/reaction_system_with_concentrations_param.py @@ -11,6 +11,9 @@ class ParameterObj(object): def __repr__(self): return "@{0}".format(self.name) +def is_param(some_object): + if isinstance(some_object, ParameterObj): + return True class ReactionSystemWithConcentrationsParam(ReactionSystem): @@ -251,9 +254,21 @@ class ReactionSystemWithConcentrationsParam(ReactionSystem): " ) 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_param_reactions(self, soft=False): + print(C_MARK_INFO + " Parametric reactions:") + if soft and len(self.parametric_reactions) > 50: + print(" -> there are more than 50 reactions (" + str(len(self.parametric_reactions)) + ")") + else: + print(" "*4 + "{0: ^35}{1: ^25}{2: ^15}".format("reactants"," inhibitors"," products")) + for reactants,inhibitors,products in self.parametric_reactions: + print(" " + "- {0: ^35}{1: ^25}{2: ^15}".format( + "{ " + self.state_to_str(reactants) + " }", + " { " + self.state_to_str(inhibitors) + " }", + " { " + self.state_to_str(products) + " }")) def show_max_concentrations(self): - print(C_MARK_INFO + " Maximal allowed concentration levels (for optimized translation to RS):") + 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)) @@ -265,6 +280,7 @@ class ReactionSystemWithConcentrationsParam(ReactionSystem): 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() diff --git a/smt/smt_checker_rsc_param.py b/smt/smt_checker_rsc_param.py index b0ec3bb..e6b2d07 100644 --- a/smt/smt_checker_rsc_param.py +++ b/smt/smt_checker_rsc_param.py @@ -11,6 +11,8 @@ from colour import * from logics import rsLTL_Encoder +from rs.reaction_system_with_concentrations_param import ParameterObj, is_param + # def simplify(x): # return x @@ -157,19 +159,22 @@ class SmtCheckerRSCParam(object): background set. """ - # TODO: - # - # We should be checking here if we acutally need - # any variables for parameters, because there might - # be no parameters in RS - level = self.next_level_to_encode - - params = [] - for entity in self.rs.background_set: - new_var = Int("L{:d}_Pm_{:s}".format(level, entity)) - params.append(new_var) - self.v_param.append(params) + + param_vars_for_cur_level = dict() + + 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.background_set: + new_var = Int("L{:d}_Pm_{:s}".format(level, entity)) + vars_for_param.append(new_var) + + param_vars_for_cur_level[param_name] = vars_for_param + + self.v_param.append(param_vars_for_cur_level) print(self.v_param) @@ -270,6 +275,21 @@ class SmtCheckerRSCParam(object): return enc_reaction + def enc_single_param_reaction(self, level, p_reaction): + + reactants, inhibitors, products = p_reaction + + if is_param(reactants): + print("PARAM") + + if is_param(inhibitors): + print("PARAM") + + if is_param(products): + print("PARAM") + + return True + def enc_rs_trans(self, level): """Encodes the transition relation""" @@ -293,6 +313,11 @@ class SmtCheckerRSCParam(object): enc_reaction = self.enc_single_reaction(level, reaction) enc_trans = simplify(And(enc_trans, enc_reaction)) + for p_reaction in self.rs.parametric_reactions: + print(p_reaction) + enc_p_reaction = self.enc_single_param_reaction(level, p_reaction) + enc_trans = simplify(And(enc_trans, enc_p_reaction)) + # Next we encode the MAX concentration values: # we collect those from the intermediate product variables From 9297a7caf7f1e33a655d240096bb501529012fcf Mon Sep 17 00:00:00 2001 From: Artur Meski Date: Sun, 17 Sep 2017 12:24:07 +0100 Subject: [PATCH 23/35] Semi-working version of parametric MC for rsLTL --- ...action_system_with_concentrations_param.py | 23 +- rs_testing.py | 38 ++- rssmt.py | 16 +- smt/smt_checker_rsc_param.py | 228 +++++++++++------- 4 files changed, 188 insertions(+), 117 deletions(-) diff --git a/rs/reaction_system_with_concentrations_param.py b/rs/reaction_system_with_concentrations_param.py index c06bdc3..38dcf11 100644 --- a/rs/reaction_system_with_concentrations_param.py +++ b/rs/reaction_system_with_concentrations_param.py @@ -14,6 +14,8 @@ class ParameterObj(object): def is_param(some_object): if isinstance(some_object, ParameterObj): return True + else: + return False class ReactionSystemWithConcentrationsParam(ReactionSystem): @@ -21,7 +23,7 @@ class ReactionSystemWithConcentrationsParam(ReactionSystem): self.reactions = [] self.parameters = dict() - self.parametric_reactions = [] + # self.parametric_reactions = [] self.meta_reactions = dict() self.permanent_entities = dict() self.background_set = [] @@ -170,10 +172,7 @@ class ReactionSystemWithConcentrationsParam(ReactionSystem): raise RuntimeError("No products defined") reaction = self.process_rip(R,I,P) - if self.is_parametric_reaction(reaction): - self.parametric_reactions.append(reaction) - else: - self.reactions.append(reaction) + self.reactions.append(reaction) def add_reaction_without_reactants(self, R, I, P): """Adds a reaction""" @@ -254,18 +253,6 @@ class ReactionSystemWithConcentrationsParam(ReactionSystem): " ) 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_param_reactions(self, soft=False): - print(C_MARK_INFO + " Parametric reactions:") - if soft and len(self.parametric_reactions) > 50: - print(" -> there are more than 50 reactions (" + str(len(self.parametric_reactions)) + ")") - else: - print(" "*4 + "{0: ^35}{1: ^25}{2: ^15}".format("reactants"," inhibitors"," products")) - for reactants,inhibitors,products in self.parametric_reactions: - print(" " + "- {0: ^35}{1: ^25}{2: ^15}".format( - "{ " + self.state_to_str(reactants) + " }", - " { " + self.state_to_str(inhibitors) + " }", - " { " + self.state_to_str(products) + " }")) def show_max_concentrations(self): print(C_MARK_INFO + " Maximal allowed concentration levels:") @@ -280,7 +267,7 @@ class ReactionSystemWithConcentrationsParam(ReactionSystem): def show(self, soft=False): self.show_background_set() self.show_reactions(soft) - self.show_param_reactions(soft) + # self.show_param_reactions(soft) self.show_permanent_entities() self.show_meta_reactions() self.show_max_concentrations() diff --git a/rs_testing.py b/rs_testing.py index 4d746c6..5e1a818 100644 --- a/rs_testing.py +++ b/rs_testing.py @@ -13,11 +13,42 @@ def run_tests(): # scalable_chain(print_system=True) # example44() # example44_param() - trivial_param() - + simple_param() def trivial_param(): + r = ReactionSystemWithConcentrationsParam() + r.add_bg_set_entity(("x", 3)) + r.add_bg_set_entity(("c", 3)) + r.add_bg_set_entity(("final", 1)) + + param_p1 = r.get_param("P1") + + r.add_reaction(param_p1, [("c", 1)], [("final", 1)]) + # r.add_reaction([("x", 1)], [("c", 1)], [("final", 1)]) + + c = ContextAutomatonWithConcentrations(r) + c.add_init_state("0") + c.add_state("1") + c.add_transition("0", [("x", 1)], "1") + c.add_transition("1", [], "1") + + rc = ReactionSystemWithAutomaton(r, c) + rc.show() + smt_rsc = SmtCheckerRSCParam(rc) + + f1 = 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(): + r = ReactionSystemWithConcentrationsParam() r.add_bg_set_entity(("x", 3)) r.add_bg_set_entity(("y", 3)) @@ -29,7 +60,8 @@ def trivial_param(): r.add_reaction([("x", 1)], [("c", 1)], [("y", 2)]) r.add_reaction([("y", 1)], [("c", 1)], [("z", 1)]) - r.add_reaction(param_p1, [("c", 1)], [("final", 1)]) + # r.add_reaction([("y", 1)], [("c", 1)], r.get_param("P2")) + r.add_reaction(param_p1, [("x", 1)], [("final", 1)]) c = ContextAutomatonWithConcentrations(r) c.add_init_state("0") diff --git a/rssmt.py b/rssmt.py index e37f75e..05cca65 100755 --- a/rssmt.py +++ b/rssmt.py @@ -18,7 +18,7 @@ import resource profiling = False -################################################################## +################################################################## version = "2017/04/08/00" rsmc_banner = """ @@ -28,24 +28,26 @@ Version: """ + version + """ Author: Artur Męski / """ -################################################################## +################################################################## + 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() -################################################################## - +################################################################## + if __name__ == "__main__": try: if profiling: diff --git a/smt/smt_checker_rsc_param.py b/smt/smt_checker_rsc_param.py index e6b2d07..3520d00 100644 --- a/smt/smt_checker_rsc_param.py +++ b/smt/smt_checker_rsc_param.py @@ -13,8 +13,8 @@ from logics import rsLTL_Encoder from rs.reaction_system_with_concentrations_param import ParameterObj, is_param -# def simplify(x): -# return x +def simplify(x): + return x def z3_max(a, b): return If(a > b, a, b) @@ -44,7 +44,8 @@ class SmtCheckerRSCParam(object): self.v_improd = [] self.v_improd_for_entities = [] - self.v_param = [] + # parameters: + self.v_param = dict() self.next_level_to_encode = 0 @@ -58,10 +59,10 @@ class SmtCheckerRSCParam(object): # improducible - entities that are never produces (there is no reaction that produces that entity) self.loop_position = Int("loop_position") - self.solver = Solver() #For("QF_FD") - self.verification_time = None + + self.prepare_param_variables() def reset(self): """Reinitialises the state of the checker""" @@ -74,7 +75,6 @@ class SmtCheckerRSCParam(object): self.prepare_state_variables() self.prepare_context_variables() self.prepare_intermediate_product_variables() - self.prepare_param_variables() self.next_level_to_encode += 1 def prepare_context_variables(self): @@ -133,22 +133,33 @@ class SmtCheckerRSCParam(object): reaction_id = self.rs.reactions.index(reaction) entities_dict = dict() - for entity, conc in products: - new_var = Int("IP" + str(level) + "_R" + - str(reaction_id) + "_e" + str(entity)) - entities_dict[entity] = new_var + + if is_param(products): - all_entities_dict.setdefault(entity, []) - all_entities_dict[entity].append(new_var) + for entity in self.rs.set_of_bgset_ids: + entity_name = self.rs.get_entity_name(entity) + new_var = Int("L"+ str(level) + "_ImProd" + "_r" + + str(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("L"+ str(level) + "_ImProd" + "_r" + + str(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.v_improd.append(reactions_dict) self.v_improd_for_entities.append(all_entities_dict) - - # print(self.v_improd) - - # print(self.v_improd_for_entities) def prepare_param_variables(self): """ @@ -159,42 +170,55 @@ class SmtCheckerRSCParam(object): background set. """ - level = self.next_level_to_encode - - param_vars_for_cur_level = dict() - 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.background_set: - new_var = Int("L{:d}_Pm_{:s}".format(level, entity)) + 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) - param_vars_for_cur_level[param_name] = vars_for_param - - self.v_param.append(param_vars_for_cur_level) - - print(self.v_param) - + self.v_param[param_name] = vars_for_param + + def enc_param_concentration_levels_assertion(self): + """ + Assertions for the parameter variables + """ + enc_param_gz = True + enc_param_at_least_one = False + for param_vars in self.v_param.values(): + for pvar in param_vars: + # + # TODO: fixed upper limit: 100 + # + 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)) + + return simplify(And(enc_param_gz, enc_param_at_least_one)) + def enc_concentration_levels_assertion(self, level): """ - Encodes assertions that (some) variables need to be >0 + 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] + enc_gz = True + + for e_i in self.rs.set_of_bgset_ids: + var = self.v[level][e_i] + var_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)) + enc_gz = simplify(And(enc_gz, var >= 0, var_ctx >= 0, var <= e_max, var_ctx <= e_max)) + + vars_per_reaction = self.v_improd_for_entities[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_nz + return enc_gz def enc_init_state(self, level): """Encodes the initial state at the given level""" @@ -233,39 +257,45 @@ class SmtCheckerRSCParam(object): # we need reaction_id to find the intermediate product variable reaction_id = self.rs.reactions.index(reaction) - # ** REACTANTS ** - + # ** REACTANTS ******************************************* enc_reactants = True - for entity, conc in reactants: - enc_reactants = And(enc_reactants, Or( - self.v[level][entity] >= conc, self.v_ctx[level][entity] >= conc)) - - # ** INHIBITORS ** + 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[level][entity] >= self.v_param[param_name][entity], + self.v_ctx[level][entity] >= self.v_param[param_name][entity])) + else: + for entity, conc in reactants: + enc_reactants = And(enc_reactants, Or( + self.v[level][entity] >= conc, + self.v_ctx[level][entity] >= conc)) + # ** INHIBITORS ****************************************** enc_inhibitors = True - for entity, conc in inhibitors: - enc_inhibitors = And(enc_inhibitors, And( - self.v[level][entity] < conc, self.v_ctx[level][entity] < conc)) + 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[level][entity] < self.v_param[param_name][entity], + self.v_ctx[level][entity] < self.v_param[param_name][entity])) + else: + for entity, conc in inhibitors: + enc_inhibitors = And(enc_inhibitors, And( + self.v[level][entity] < conc, + self.v_ctx[level][entity] < conc)) - # ** PRODUCTS ** - - # entities that were produced - produced = [] + # ** PRODUCTS ******************************************* enc_products = True - for entity, conc in products: - enc_products = And(enc_products, self.v_improd[ - level + 1][reaction_id][entity] == conc) - produced.append(entity) - - print("----", str(produced)) - - # encoding of the zeros (for the products) - # we iterate through all the entities' indices - - # for entity in range(len(self.rs.background_set)): - # if entity not in produced: - # enc_products = And(enc_products, self.v_improd[ - # level + 1][reaction_id][entity] == 0) + if is_param(products): + param_name = products.name + for entity in self.rs.set_of_bgset_ids: + enc_products = And(enc_products, + self.v_improd[level + 1][reaction_id][entity] == self.v_param[param_name][entity]) + else: + for entity, conc in products: + enc_products = And(enc_products, + self.v_improd[level + 1][reaction_id][entity] == conc) # # (R and I) iff P @@ -275,20 +305,25 @@ class SmtCheckerRSCParam(object): return enc_reaction - def enc_single_param_reaction(self, level, p_reaction): + def enc_general_reaction_enabledness(self, level): + """ + General enabledness condition for reactions - reactants, inhibitors, products = p_reaction + 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. - if is_param(reactants): - print("PARAM") + 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) + """ - if is_param(inhibitors): - print("PARAM") - - if is_param(products): - print("PARAM") - - return True + enc_cond = False + for entity in self.rs.set_of_bgset_ids: + enc_cond = simplify(Or(enc_cond, self.v[level][entity] > 0, self.v_ctx[level][entity] > 0)) + + return enc_cond def enc_rs_trans(self, level): """Encodes the transition relation""" @@ -313,11 +348,6 @@ class SmtCheckerRSCParam(object): enc_reaction = self.enc_single_reaction(level, reaction) enc_trans = simplify(And(enc_trans, enc_reaction)) - for p_reaction in self.rs.parametric_reactions: - print(p_reaction) - enc_p_reaction = self.enc_single_param_reaction(level, p_reaction) - enc_trans = simplify(And(enc_trans, enc_p_reaction)) - # Next we encode the MAX concentration values: # we collect those from the intermediate product variables @@ -334,17 +364,16 @@ class SmtCheckerRSCParam(object): # {products}[level+1] # current_v_improd_for_entities = self.v_improd_for_entities[level + 1] - for entity, per_reaction_vars in current_v_improd_for_entities.items(): - enc_max_prod = simplify( And(enc_max_prod, self.v[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) - # for entity in self.improducible_entities: - # enc_max_prod = simplify( - # And(enc_max_prod, self.v[level + 1][entity] == 0)) + enc_trans_with_max = simplify(And(enc_general_cond, enc_max_prod, enc_trans)) - enc_trans_with_max = simplify(And(enc_max_prod, enc_trans)) + print(enc_trans_with_max) return enc_trans_with_max @@ -482,10 +511,30 @@ class SmtCheckerRSCParam(object): " " + 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(" }") + def check_rsltl( - self, formula, print_witness=True, print_time=True, print_mem=True, + self, formula, + print_witness=True, print_time=True, print_mem=True, max_level=None): """Bounded Model Checking for rsLTL properties""" @@ -506,6 +555,7 @@ class SmtCheckerRSCParam(object): self.prepare_all_variables() self.solver.add(self.enc_concentration_levels_assertion(0)) + self.solver.add(self.enc_param_concentration_levels_assertion()) encoder = rsLTL_Encoder(self) From a43277c159e63d39af7333da3a9874ec083462d0 Mon Sep 17 00:00:00 2001 From: Artur Meski Date: Sun, 17 Sep 2017 17:00:21 +0100 Subject: [PATCH 24/35] Fixing the encoding of reactions --- smt/smt_checker_rsc_param.py | 72 +++++++++++++++++++----------------- 1 file changed, 39 insertions(+), 33 deletions(-) diff --git a/smt/smt_checker_rsc_param.py b/smt/smt_checker_rsc_param.py index 3520d00..b2adc14 100644 --- a/smt/smt_checker_rsc_param.py +++ b/smt/smt_checker_rsc_param.py @@ -13,8 +13,8 @@ from logics import rsLTL_Encoder from rs.reaction_system_with_concentrations_param import ParameterObj, is_param -def simplify(x): - return x +# def simplify(x): +# return x def z3_max(a, b): return If(a > b, a, b) @@ -263,6 +263,7 @@ class SmtCheckerRSCParam(object): 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.v[level][entity] >= self.v_param[param_name][entity], self.v_ctx[level][entity] >= self.v_param[param_name][entity])) else: @@ -277,8 +278,10 @@ class SmtCheckerRSCParam(object): param_name = inhibitors.name for entity in self.rs.set_of_bgset_ids: enc_inhibitors = And(enc_inhibitors, Or( - self.v[level][entity] < self.v_param[param_name][entity], - self.v_ctx[level][entity] < self.v_param[param_name][entity])) + self.v_param[param_name][entity] == 0, + And( + self.v[level][entity] < self.v_param[param_name][entity], + self.v_ctx[level][entity] < self.v_param[param_name][entity]))) else: for entity, conc in inhibitors: enc_inhibitors = And(enc_inhibitors, And( @@ -290,18 +293,35 @@ class SmtCheckerRSCParam(object): if is_param(products): param_name = products.name for entity in self.rs.set_of_bgset_ids: - enc_products = And(enc_products, - self.v_improd[level + 1][reaction_id][entity] == self.v_param[param_name][entity]) + enc_products = simplify(And(enc_products, Or( + self.v_param[param_name][entity] == 0, + self.v_improd[level + 1][reaction_id][entity] == self.v_param[param_name][entity]))) else: for entity, conc in products: - enc_products = And(enc_products, - self.v_improd[level + 1][reaction_id][entity] == conc) + enc_products = simplify(And(enc_products, + self.v_improd[level + 1][reaction_id][entity] == conc)) + enc_no_prod = True + if is_param(products): + for entity in self.rs.set_of_bgset_ids: + enc_products = And(enc_products, + self.v_improd[level + 1][reaction_id][entity] == 0) + else: + for entity, _ in products: + enc_no_prod = simplify(And(enc_no_prod, + self.v_improd[level + 1][reaction_id][entity] == 0)) + # # (R and I) iff P # - enc_reaction = simplify( - And(enc_reactants, enc_inhibitors) == enc_products) + 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 @@ -364,7 +384,8 @@ class SmtCheckerRSCParam(object): # {products}[level+1] # current_v_improd_for_entities = self.v_improd_for_entities[level + 1] - for entity, per_reaction_vars in current_v_improd_for_entities.items(): + 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.v[level + 1][entity] == self.enc_max(per_reaction_vars))) @@ -373,15 +394,18 @@ class SmtCheckerRSCParam(object): enc_trans_with_max = simplify(And(enc_general_cond, enc_max_prod, enc_trans)) - print(enc_trans_with_max) + # print(enc_trans_with_max) return enc_trans_with_max def enc_max(self, elements): enc = None + + if elements == []: + enc = 0 - if len(elements) == 1: + elif len(elements) == 1: enc = z3_max(0, elements[0]) elif len(elements) > 1: @@ -425,22 +449,6 @@ class SmtCheckerRSCParam(object): raise RuntimeError("Should not be used with RSC") - # enc = True - # used_entities_ids = self.rs.get_state_ids(state) - # - # for ent,conc in state: - # e_id = self.rs.get_entity_id(ent) - # enc = And(enc, self.v[level][e_id] == conc) - # - # not_in_state = set(range(len(self.rs.background_set))) - # not_in_state = not_in_state.difference(set(used_entities_ids)) - # - # for entity in not_in_state: - # enc = And(enc, self.v[level][entity] == 0) - # - # return simplify(enc) - - def enc_min_state(self, level, state): """Encodes the state at the given level with the minimal required concentration levels""" @@ -454,8 +462,7 @@ class SmtCheckerRSCParam(object): # for entity in state_ids: # enc = And(enc, self.v[level][entity]) - return simplify(enc) - + return simplify(enc) def enc_state_with_blocking(self, level, prop): """Encodes the state at the given level with blocking certain concentrations""" @@ -472,7 +479,6 @@ class SmtCheckerRSCParam(object): enc = And(enc, self.v[level][e_id] < conc) return simplify(enc) - def decode_witness(self, max_level, print_model=False): @@ -597,7 +603,7 @@ class SmtCheckerRSCParam(object): "[" + colour_str(C_BOLD, "+") + "] " + colour_str( C_GREEN, "SAT at level=" + str(self.current_level))) - print(self.solver.model()) + # print(self.solver.model()) if print_witness: print("\n{:=^70}".format("[ WITNESS ]")) self.decode_witness(self.current_level) From a3f06724bf98883c64c02cd66d12604b1c87c0d6 Mon Sep 17 00:00:00 2001 From: Artur Meski Date: Sun, 17 Sep 2017 19:02:46 +0100 Subject: [PATCH 25/35] Production for parameters --- smt/smt_checker_rsc_param.py | 61 ++++++++++++++++++++++++++---------- 1 file changed, 45 insertions(+), 16 deletions(-) diff --git a/smt/smt_checker_rsc_param.py b/smt/smt_checker_rsc_param.py index b2adc14..02b8231 100644 --- a/smt/smt_checker_rsc_param.py +++ b/smt/smt_checker_rsc_param.py @@ -60,6 +60,7 @@ class SmtCheckerRSCParam(object): self.loop_position = Int("loop_position") self.solver = Solver() #For("QF_FD") + self.verification_time = None self.prepare_param_variables() @@ -185,18 +186,26 @@ class SmtCheckerRSCParam(object): """ Assertions for the parameter variables """ + + if len(self.v_param) == 0: + return True + enc_param_gz = True - enc_param_at_least_one = False + 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 - # + + # 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)) - - return simplify(And(enc_param_gz, enc_param_at_least_one)) - + + enc_non_empty = simplify(And(enc_non_empty, enc_param_at_least_one)) + + return simplify(And(enc_param_gz, enc_non_empty)) + def enc_concentration_levels_assertion(self, level): """ Encodes assertions that (some) variables need to be >=0 @@ -293,18 +302,18 @@ class SmtCheckerRSCParam(object): if is_param(products): param_name = products.name for entity in self.rs.set_of_bgset_ids: - enc_products = simplify(And(enc_products, Or( - self.v_param[param_name][entity] == 0, - self.v_improd[level + 1][reaction_id][entity] == self.v_param[param_name][entity]))) + enc_products = simplify(And(enc_products, + self.v_improd[level + 1][reaction_id][entity] == self.v_param[param_name][entity])) else: for entity, conc in products: enc_products = simplify(And(enc_products, self.v_improd[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_products = And(enc_products, + enc_no_prod = And(enc_no_prod, self.v_improd[level + 1][reaction_id][entity] == 0) else: for entity, _ in products: @@ -481,6 +490,11 @@ class SmtCheckerRSCParam(object): return simplify(enc) def decode_witness(self, max_level, print_model=False): + """ + Decodes the witness + + Also decodes the parameters + """ m = self.solver.model() @@ -537,12 +551,26 @@ class SmtCheckerRSCParam(object): " " + str(self.rs.get_entity_name(entity)) + "=" + str(var_rep), end="") print(" }") + + print("\n") def check_rsltl( self, formula, - print_witness=True, print_time=True, print_mem=True, - max_level=None): - """Bounded Model Checking for rsLTL properties""" + 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) + """ self.reset() @@ -603,11 +631,12 @@ class SmtCheckerRSCParam(object): "[" + colour_str(C_BOLD, "+") + "] " + colour_str( C_GREEN, "SAT at level=" + str(self.current_level))) - # print(self.solver.model()) + print(self.solver.model()) if print_witness: print("\n{:=^70}".format("[ WITNESS ]")) self.decode_witness(self.current_level) - break + if not cont_if_sat: + break else: self.solver.pop() From 9501feeb7ed33fa3679a0318b1850eda27d97e0e Mon Sep 17 00:00:00 2001 From: Artur Meski Date: Sun, 17 Sep 2017 19:53:08 +0100 Subject: [PATCH 26/35] Testing example --- rs_testing.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/rs_testing.py b/rs_testing.py index 5e1a818..f488eac 100644 --- a/rs_testing.py +++ b/rs_testing.py @@ -59,9 +59,10 @@ def simple_param(): 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([("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") @@ -80,7 +81,7 @@ def simple_param(): # # WARNING: depth limit is set # - smt_rsc.check_rsltl(formula=f1, max_level=10, print_witness=True) + smt_rsc.check_rsltl(formula=f1, max_level=10, print_witness=True, cont_if_sat=False) def example44_param(): From 687b1c9cab58ab6d4024b5eb23903206752ce230 Mon Sep 17 00:00:00 2001 From: Artur Meski Date: Sun, 17 Sep 2017 20:57:45 +0100 Subject: [PATCH 27/35] maxsat --- rs_testing.py | 2 +- smt/smt_checker_rsc_param.py | 64 +++++++++++++++++++++++------------- 2 files changed, 42 insertions(+), 24 deletions(-) diff --git a/rs_testing.py b/rs_testing.py index f488eac..8511320 100644 --- a/rs_testing.py +++ b/rs_testing.py @@ -72,7 +72,7 @@ def simple_param(): rc = ReactionSystemWithAutomaton(r, c) rc.show() - smt_rsc = SmtCheckerRSCParam(rc) + smt_rsc = SmtCheckerRSCParam(rc, optimise=True) f1 = Formula_rsLTL.f_F( BagDescription.f_TRUE(), diff --git a/smt/smt_checker_rsc_param.py b/smt/smt_checker_rsc_param.py index 02b8231..77a3683 100644 --- a/smt/smt_checker_rsc_param.py +++ b/smt/smt_checker_rsc_param.py @@ -21,7 +21,7 @@ def z3_max(a, b): class SmtCheckerRSCParam(object): - def __init__(self, rsca): + def __init__(self, rsca, optimise=False): rsca.sanity_check() @@ -30,6 +30,8 @@ class SmtCheckerRSCParam(object): self.rs = rsca.rs self.ca = rsca.ca + + self.optimise = optimise self.initialise() @@ -59,7 +61,11 @@ class SmtCheckerRSCParam(object): # improducible - entities that are never produces (there is no reaction that produces that entity) self.loop_position = Int("loop_position") - self.solver = Solver() #For("QF_FD") + + if self.optimise: + self.solver = Optimize() + else: + self.solver = Solver() #For("QF_FD") self.verification_time = None @@ -205,6 +211,12 @@ class SmtCheckerRSCParam(object): 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): """ @@ -452,7 +464,6 @@ class SmtCheckerRSCParam(object): return enc_trans - def enc_exact_state(self, level, state): """Encodes the state at the given level with the exact concentration values""" @@ -583,19 +594,22 @@ class SmtCheckerRSCParam(object): start = resource.getrusage(resource.RUSAGE_SELF).ru_utime self.prepare_all_variables() - self.solver.add(self.enc_init_state(0)) + 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)) - self.solver.add(self.enc_param_concentration_levels_assertion()) + self.solver_add(self.enc_concentration_levels_assertion(0)) + self.solver_add(self.enc_param_concentration_levels_assertion()) + + if self.optimise: + self.assert_param_optimisation() encoder = rsLTL_Encoder(self) while True: self.prepare_all_variables() - self.solver.add( + self.solver_add( self.enc_concentration_levels_assertion( self.current_level + 1)) @@ -619,11 +633,11 @@ class SmtCheckerRSCParam(object): "] Adding the formula to the solver...") encoder.flush_cache() - self.solver.add(f) + self.solver_add(f) print("[" + colour_str(C_BOLD, "i") + "] Adding the loops encoding...") - self.solver.add(self.get_loop_encodings()) + self.solver_add(self.get_loop_encodings()) result = self.solver.check() if result == sat: @@ -642,7 +656,7 @@ class SmtCheckerRSCParam(object): print("[" + colour_str(C_BOLD, "i") + "] Unrolling the transition relation") - self.solver.add(self.enc_transition_relation(self.current_level)) + self.solver_add(self.enc_transition_relation(self.current_level)) print( "{:->70}".format("[ level=" + str(self.current_level) + " done ]")) @@ -668,7 +682,6 @@ class SmtCheckerRSCParam(object): repr( resource.getrusage(resource.RUSAGE_SELF).ru_maxrss / (1024 * 1024)) + " MB")) - def dummy_unroll(self, levels): """Unrolls the variables for testing purposes""" @@ -680,7 +693,6 @@ class SmtCheckerRSCParam(object): print(C_MARK_INFO + " Dummy Unrolling done.") - def state_equality(self, level_A, level_B): """Encodes equality of two states at two different levels""" @@ -693,8 +705,7 @@ class SmtCheckerRSCParam(object): 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 - + return eq_enc def get_loop_encodings(self): @@ -713,8 +724,15 @@ class SmtCheckerRSCParam(object): loop_enc = simplify(And(loop_enc, Implies( loop_var == i, self.state_equality(i-1, k) ))) return loop_enc + + def solver_add(self, expression): + + if expression == False: + raise RuntimeError("Trying to assert False.") - + if not (expression == True): + self.solver.add(expression) + def check_reachability(self, state, print_witness=True, print_time=True, print_mem=True, max_level=1000): """Main testing function""" @@ -726,16 +744,16 @@ class SmtCheckerRSCParam(object): start = resource.getrusage(resource.RUSAGE_SELF).ru_utime self.prepare_all_variables() - self.solver.add(self.enc_init_state(0)) + 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)) + 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)) + 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() @@ -744,7 +762,7 @@ class SmtCheckerRSCParam(object): 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)) + self.solver_add(self.enc_state_with_blocking(self.current_level,state)) result = self.solver.check() if result == sat: @@ -757,7 +775,7 @@ class SmtCheckerRSCParam(object): self.solver.pop() print("[" + colour_str(C_BOLD, "i") + "] Unrolling the transition relation") - self.solver.add(self.enc_transition_relation(self.current_level)) + self.solver_add(self.enc_transition_relation(self.current_level)) print("{:->70}".format("[ level=" + str(self.current_level) + " done ]")) self.current_level += 1 @@ -790,7 +808,7 @@ class SmtCheckerRSCParam(object): self.prepare_all_variables() init_s = self.enc_init_state(0) print(init_s) - self.solver.add(init_s) + self.solver_add(init_s) self.current_level = 0 self.prepare_all_variables() @@ -808,7 +826,7 @@ class SmtCheckerRSCParam(object): s = self.enc_min_state(self.current_level,state) print("Test: ", s) - self.solver.add(s) + self.solver_add(s) result = self.solver.check() if result == sat: @@ -822,7 +840,7 @@ class SmtCheckerRSCParam(object): print("[i] Unrolling the transition relation") t = self.enc_transition_relation(self.current_level) print(t) - self.solver.add(t) + self.solver_add(t) print("-----[ level=" + str(self.current_level) + " done ]") self.current_level += 1 From 81b1ab45076780a5e453fa47464ef89412fe3cae Mon Sep 17 00:00:00 2001 From: Artur Meski Date: Fri, 22 Sep 2017 21:22:38 +0200 Subject: [PATCH 28/35] Command line arguments --- rs_testing.py | 9 +++++---- rssmt.py | 18 +++++++++++++++--- smt/smt_checker_rsc_param.py | 9 +++++++-- 3 files changed, 27 insertions(+), 9 deletions(-) diff --git a/rs_testing.py b/rs_testing.py index 8511320..b0ac971 100644 --- a/rs_testing.py +++ b/rs_testing.py @@ -5,7 +5,7 @@ from logics import * import sys import resource -def run_tests(): +def run_tests(cmd_args): # test_extended_automaton() # process() @@ -13,7 +13,7 @@ def run_tests(): # scalable_chain(print_system=True) # example44() # example44_param() - simple_param() + simple_param(cmd_args) def trivial_param(): @@ -47,7 +47,7 @@ def trivial_param(): smt_rsc.check_rsltl(formula=f1, max_level=10, print_witness=True) -def simple_param(): +def simple_param(cmd_args): r = ReactionSystemWithConcentrationsParam() r.add_bg_set_entity(("x", 3)) @@ -72,7 +72,8 @@ def simple_param(): rc = ReactionSystemWithAutomaton(r, c) rc.show() - smt_rsc = SmtCheckerRSCParam(rc, optimise=True) + + smt_rsc = SmtCheckerRSCParam(rc, optimise=cmd_args.optimise) f1 = Formula_rsLTL.f_F( BagDescription.f_TRUE(), diff --git a/rssmt.py b/rssmt.py index 05cca65..c3a1126 100755 --- a/rssmt.py +++ b/rssmt.py @@ -6,6 +6,8 @@ """ +import argparse + from rs import * from smt import * import sys @@ -14,10 +16,11 @@ import rs_testing from colour import * -import resource - profiling = False +if profiling: + import resource + ################################################################## version = "2017/04/08/00" @@ -43,8 +46,17 @@ def print_banner(): def main(): """Main function""" + parser = argparse.ArgumentParser() + + parser.add_argument("-v", "--verbose", + help="turn verbosity on", action="store_true") + parser.add_argument("-o", "--optimise", + help="minimise the parametric computation result", action="store_true") + + args = parser.parse_args() + print_banner() - rs_testing.run_tests() + rs_testing.run_tests(args) ################################################################## diff --git a/smt/smt_checker_rsc_param.py b/smt/smt_checker_rsc_param.py index 77a3683..2b6a005 100644 --- a/smt/smt_checker_rsc_param.py +++ b/smt/smt_checker_rsc_param.py @@ -726,12 +726,17 @@ class SmtCheckerRSCParam(object): 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.") - if not (expression == True): - self.solver.add(expression) + self.solver.add(expression) def check_reachability(self, state, print_witness=True, print_time=True, print_mem=True, max_level=1000): From f7097e37d7a5de6c6369fae4013d94f1699fc06a Mon Sep 17 00:00:00 2001 From: Artur Meski Date: Sun, 5 Nov 2017 17:51:04 +0000 Subject: [PATCH 29/35] Colour prints --- colour.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/colour.py b/colour.py index 29ba737..80b7c65 100644 --- a/colour.py +++ b/colour.py @@ -13,7 +13,16 @@ C_MARK_ERROR = C_BOLD + "[" + C_RED + "!" + C_ENDC + C_BOLD + "]" + C_ENDC def colour_str(col, s): return col + s + C_ENDC +def green_str(s): + return C_GREEN + s + C_ENDC + def print_error(s): print(C_MARK_ERROR + " " + C_RED + s + C_ENDC) - + +def print_info(s): + print("[" + colour_str(C_BOLD, "i") + "] {:s}".format(s)) + +def print_positive(s): + print("[" + colour_str(C_BOLD, "+") + "] {:s}".format(s)) + # EOF From d910de17ab6f3e5558e2d94c14fb141d57cdca4c Mon Sep 17 00:00:00 2001 From: Artur Meski Date: Sun, 5 Nov 2017 17:51:18 +0000 Subject: [PATCH 30/35] New tests --- rs_testing.py | 92 ++++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 91 insertions(+), 1 deletion(-) diff --git a/rs_testing.py b/rs_testing.py index b0ac971..5733617 100644 --- a/rs_testing.py +++ b/rs_testing.py @@ -13,7 +13,8 @@ def run_tests(cmd_args): # scalable_chain(print_system=True) # example44() # example44_param() - simple_param(cmd_args) + heat_shock_response_param(cmd_args) + # simple_param(cmd_args) def trivial_param(): @@ -85,6 +86,95 @@ def simple_param(cmd_args): 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_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() From 1b94cbc7c00b26327ca806cc6fd4a7db12b0eafb Mon Sep 17 00:00:00 2001 From: Artur Meski Date: Sun, 5 Nov 2017 17:52:02 +0000 Subject: [PATCH 31/35] Loading of the variables to use for the encoding --- logics/rsltl_encoder.py | 24 +++++++++++++++++++++--- 1 file changed, 21 insertions(+), 3 deletions(-) diff --git a/logics/rsltl_encoder.py b/logics/rsltl_encoder.py index 79844b8..c1f61e7 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): From a9a37baba7d6ec343fc5323e3e5756190d90b114 Mon Sep 17 00:00:00 2001 From: Artur Meski Date: Sun, 5 Nov 2017 17:52:51 +0000 Subject: [PATCH 32/35] Name --- rssmt.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rssmt.py b/rssmt.py index c3a1126..6e50d46 100755 --- a/rssmt.py +++ b/rssmt.py @@ -28,7 +28,7 @@ rsmc_banner = """ Reaction Systems SMT-Based Model Checking Version: """ + version + """ -Author: Artur Męski / +Author: Artur Meski / """ ################################################################## From 4dbeac742dca2c7ad1bcf167f6876b9f98916ad5 Mon Sep 17 00:00:00 2001 From: Artur Meski Date: Sun, 5 Nov 2017 17:53:08 +0000 Subject: [PATCH 33/35] Testing of multiple rsLTL formulae at once --- smt/smt_checker_rsc_param.py | 316 +++++++++++++++++++++++------------ 1 file changed, 207 insertions(+), 109 deletions(-) diff --git a/smt/smt_checker_rsc_param.py b/smt/smt_checker_rsc_param.py index 2b6a005..fb264ed 100644 --- a/smt/smt_checker_rsc_param.py +++ b/smt/smt_checker_rsc_param.py @@ -37,14 +37,24 @@ class SmtCheckerRSCParam(object): def initialise(self): """Initialises all the variables used by the checker""" - - self.v = [] - self.v_ctx = [] - self.ca_state = [] + + ### "Currently" used variables (loaded from self.path_v...) + self.v = None + self.v_ctx = None + self.ca_state = None # intermediate products: - self.v_improd = [] - self.v_improd_for_entities = [] + 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() @@ -60,9 +70,11 @@ class SmtCheckerRSCParam(object): # improd - intermediate products # improducible - entities that are never produces (there is no reaction that produces that entity) - self.loop_position = Int("loop_position") + # TODO: number of loops == number of paths + self.loop_position = None + self.path_loop_position = dict() - if self.optimise: + if self.optimise: self.solver = Optimize() else: self.solver = Solver() #For("QF_FD") @@ -74,40 +86,67 @@ class SmtCheckerRSCParam(object): def reset(self): """Reinitialises the state of the checker""" - self.initialise() + self.initialise() - def prepare_all_variables(self): - """Prepares all the variables""" + def prepare_all_variables(self, num_of_paths): - self.prepare_state_variables() - self.prepare_context_variables() - self.prepare_intermediate_product_variables() - self.next_level_to_encode += 1 + for path_idx in range(num_of_paths): + self.prepare_all_path_variables(path_idx) + self.next_level_to_encode += 1 - def prepare_context_variables(self): + 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: - variables.append(Int("C"+str(level)+"_"+entity)) + new_var = Int("p{:d}C{:d}_{:s}".format(path_idx, level, entity)) + variables.append(new_var) - self.v_ctx.append(variables) + self.path_v_ctx[path_idx].append(variables) - def prepare_state_variables(self): + 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: - variables.append(Int("L"+str(level)+"_"+entity)) - self.v.append(variables) + new_var = Int("p{:d}L{:d}_{:s}".format(path_idx, level, entity)) + variables.append(new_var) + self.path_v[path_idx].append(variables) - self.ca_state.append(Int("CA"+str(level)+"_state")) + # 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): + def prepare_intermediate_product_variables(self, path_idx): """ Prepares the intermediate product variables carrying the individual concentration levels produced @@ -119,14 +158,20 @@ class SmtCheckerRSCParam(object): 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.v_improd.append(None) - self.v_improd_for_entities.append(None) + self.path_v_improd[path_idx].append(None) + self.path_v_improd_for_entities[path_idx].append(None) else: reactions_dict = dict() @@ -145,8 +190,7 @@ class SmtCheckerRSCParam(object): for entity in self.rs.set_of_bgset_ids: entity_name = self.rs.get_entity_name(entity) - new_var = Int("L"+ str(level) + "_ImProd" + "_r" + - str(reaction_id) + "_" + entity_name) + 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, []) @@ -156,8 +200,7 @@ class SmtCheckerRSCParam(object): for entity, conc in products: entity_name = self.rs.get_entity_name(entity) - new_var = Int("L"+ str(level) + "_ImProd" + "_r" + - str(reaction_id) + "_" + entity_name) + 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, []) @@ -165,8 +208,8 @@ class SmtCheckerRSCParam(object): reactions_dict[reaction_id] = entities_dict - self.v_improd.append(reactions_dict) - self.v_improd_for_entities.append(all_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): """ @@ -187,6 +230,17 @@ class SmtCheckerRSCParam(object): 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): """ @@ -218,50 +272,52 @@ class SmtCheckerRSCParam(object): for pvar in param_vars: self.solver.add_soft(pvar < 1) - def enc_concentration_levels_assertion(self, level): + 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.v[level][e_i] - var_ctx = self.v_ctx[level][e_i] + 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.v_improd_for_entities[level + 1] + 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): + 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.v[level]: + 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.ca_state[level] == self.ca.get_init_state_id() + 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): + def enc_transition_relation(self, level, path_idx): return simplify( - And(self.enc_rs_trans(level), - self.enc_automaton_trans(level))) + And(self.enc_rs_trans(level, path_idx), + self.enc_automaton_trans(level, path_idx))) - def enc_single_reaction(self, level, reaction): + def enc_single_reaction(self, level, path_idx, reaction): """ Encodes a single reaction @@ -285,13 +341,13 @@ class SmtCheckerRSCParam(object): for entity in self.rs.set_of_bgset_ids: enc_reactants = And(enc_reactants, Or( self.v_param[param_name][entity] == 0, - self.v[level][entity] >= self.v_param[param_name][entity], - self.v_ctx[level][entity] >= self.v_param[param_name][entity])) + 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.v[level][entity] >= conc, - self.v_ctx[level][entity] >= conc)) + self.path_v[path_idx][level][entity] >= conc, + self.path_v_ctx[path_idx][level][entity] >= conc)) # ** INHIBITORS ****************************************** enc_inhibitors = True @@ -301,13 +357,13 @@ class SmtCheckerRSCParam(object): enc_inhibitors = And(enc_inhibitors, Or( self.v_param[param_name][entity] == 0, And( - self.v[level][entity] < self.v_param[param_name][entity], - self.v_ctx[level][entity] < self.v_param[param_name][entity]))) + 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.v[level][entity] < conc, - self.v_ctx[level][entity] < conc)) + self.path_v[path_idx][level][entity] < conc, + self.path_v_ctx[path_idx][level][entity] < conc)) # ** PRODUCTS ******************************************* enc_products = True @@ -315,22 +371,22 @@ class SmtCheckerRSCParam(object): param_name = products.name for entity in self.rs.set_of_bgset_ids: enc_products = simplify(And(enc_products, - self.v_improd[level + 1][reaction_id][entity] == self.v_param[param_name][entity])) + 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.v_improd[level + 1][reaction_id][entity] == conc)) + 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.v_improd[level + 1][reaction_id][entity] == 0) + 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.v_improd[level + 1][reaction_id][entity] == 0)) + self.path_v_improd[path_idx][level + 1][reaction_id][entity] == 0)) # # (R and I) iff P @@ -346,7 +402,7 @@ class SmtCheckerRSCParam(object): return enc_reaction - def enc_general_reaction_enabledness(self, level): + def enc_general_reaction_enabledness(self, level, path_idx): """ General enabledness condition for reactions @@ -362,11 +418,11 @@ class SmtCheckerRSCParam(object): enc_cond = False for entity in self.rs.set_of_bgset_ids: - enc_cond = simplify(Or(enc_cond, self.v[level][entity] > 0, self.v_ctx[level][entity] > 0)) + 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): + def enc_rs_trans(self, level, path_idx): """Encodes the transition relation""" # @@ -386,7 +442,7 @@ class SmtCheckerRSCParam(object): enc_trans = True for reaction in self.rs.reactions: - enc_reaction = self.enc_single_reaction(level, reaction) + 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: @@ -404,14 +460,14 @@ class SmtCheckerRSCParam(object): # => # {products}[level+1] # - current_v_improd_for_entities = self.v_improd_for_entities[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.v[level + 1][entity] == self.enc_max(per_reaction_vars))) + 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) + 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)) @@ -437,14 +493,14 @@ class SmtCheckerRSCParam(object): return enc - def enc_automaton_trans(self, level): + 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.ca_state[level] == src - dst_enc = self.ca_state[level+1] == dst + 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))) @@ -454,10 +510,10 @@ class SmtCheckerRSCParam(object): ctx_enc = True for e,c in ctx: - ctx_enc = simplify(And(ctx_enc, self.v_ctx[level][e] == c)) + 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.v_ctx[level][e] == 0)) + 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)) @@ -500,7 +556,7 @@ class SmtCheckerRSCParam(object): return simplify(enc) - def decode_witness(self, max_level, print_model=False): + def decode_witness(self, max_level, path_idx, print_model=False): """ Decodes the witness @@ -517,8 +573,8 @@ class SmtCheckerRSCParam(object): 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]]) + 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") @@ -532,8 +588,8 @@ class SmtCheckerRSCParam(object): 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]]) + 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") @@ -566,7 +622,7 @@ class SmtCheckerRSCParam(object): print("\n") def check_rsltl( - self, formula, + self, formulae_list, print_witness=True, print_time=True, print_mem=True, max_level=None, cont_if_sat=False): @@ -582,81 +638,123 @@ class SmtCheckerRSCParam(object): 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() - print("[" + colour_str(C_BOLD, "i") + - "] Running rsLTL bounded model checking") - print("[" + colour_str(C_BOLD, "i") + "] Formula: " + str(formula)) + 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 = time() start = resource.getrusage(resource.RUSAGE_SELF).ru_utime - self.prepare_all_variables() - self.solver_add(self.enc_init_state(0)) + 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() + self.prepare_all_variables(num_of_paths) - self.solver_add(self.enc_concentration_levels_assertion(0)) - self.solver_add(self.enc_param_concentration_levels_assertion()) + # 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: - 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() + # stdout.flush() # reachability test: self.solver.push() - print("[" + colour_str(C_BOLD, "i") + - "] Generating the formula encoding...") + 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))) - f = encoder.get_encoding(formula, self.current_level) - ncalls = encoder.get_ncalls() + 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("[" + 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...") + print_info("Cache hits: {:d}, encode calls: {:d} (approx: {:d})".format( + encoder.get_cache_hits(), ncalls[0], ncalls[1])) - encoder.flush_cache() - self.solver_add(f) + encoder.flush_cache() + + print_info("Adding the formulae to the solver...") + + # print (enc_form) + + self.solver_add(enc_form) - print("[" + colour_str(C_BOLD, "i") + - "] Adding the loops encoding...") + 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( - "[" + colour_str(C_BOLD, "+") + "] " + - colour_str( - C_GREEN, "SAT at level=" + str(self.current_level))) - print(self.solver.model()) + print_positive(green_str( + "SAT at level={:d}".format(self.current_level))) + # print(self.solver.model()) if print_witness: - print("\n{:=^70}".format("[ WITNESS ]")) - self.decode_witness(self.current_level) + 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() - print("[" + colour_str(C_BOLD, "i") + - "] Unrolling the transition relation") - self.solver_add(self.enc_transition_relation(self.current_level)) + 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 ]")) From 2c399dc5c175347c4d5caea9e2780d319e86941e Mon Sep 17 00:00:00 2001 From: Artur Meski Date: Sun, 26 Nov 2017 20:06:58 +0000 Subject: [PATCH 34/35] RSLTL shortcuts and gene expr example --- rs/reaction_system_with_concentrations.py | 3 +- rs_testing.py | 65 +++++++++++++++++++++-- rsltl_shortcuts.py | 26 +++++++++ 3 files changed, 87 insertions(+), 7 deletions(-) create mode 100644 rsltl_shortcuts.py 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_testing.py b/rs_testing.py index 5733617..05586bc 100644 --- a/rs_testing.py +++ b/rs_testing.py @@ -2,6 +2,8 @@ from rs import * from smt import * import rs_examples from logics import * +from rsltl_shortcuts import * + import sys import resource @@ -13,9 +15,58 @@ def run_tests(cmd_args): # scalable_chain(print_system=True) # example44() # example44_param() - heat_shock_response_param(cmd_args) + # heat_shock_response_param(cmd_args) # simple_param(cmd_args) + gene_expression(cmd_args) + +def gene_expression(cmd_args): + """ + Simple gene expression example + """ + r = ReactionSystemWithConcentrationsParam() + + r.add_bg_set_entity(("x", 1)) + r.add_bg_set_entity(("xp", 1)) + r.add_bg_set_entity(("X", 1)) + + r.add_bg_set_entity(("y", 1)) + r.add_bg_set_entity(("yp", 1)) + r.add_bg_set_entity(("Y", 1)) + + r.add_bg_set_entity(("z", 1)) + r.add_bg_set_entity(("zp", 1)) + r.add_bg_set_entity(("Z", 1)) + + r.add_bg_set_entity(("h", 1)) + r.add_bg_set_entity(("Q", 1)) + r.add_bg_set_entity(("U", 1)) + + r.add_reaction([("x",1)],[("h",1)],[("x",1)]) + r.add_reaction([("x",1)],[("h",1)],[("xp",1)]) + r.add_reaction([("x",1),("xp",1)],[("h",1)],[("X",1)]) + 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() + + def trivial_param(): r = ReactionSystemWithConcentrationsParam() @@ -38,9 +89,11 @@ def trivial_param(): rc.show() smt_rsc = SmtCheckerRSCParam(rc) - f1 = Formula_rsLTL.f_F( - BagDescription.f_TRUE(), - BagDescription.f_entity("final") >= 1) + f1 = ltl_F(bag_entity("final") >= 1) + + # f1 = Formula_rsLTL.f_F( + # BagDescription.f_TRUE(), + # BagDescription.f_entity("final") >= 1) # # WARNING: depth limit is set @@ -152,7 +205,9 @@ def heat_shock_response_param(cmd_args, print_system=True): # 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 = Formula_rsLTL.f_F(BagDescription.f_TRUE(), (BagDescription.f_entity("mfp") > 0) ) + f_reach_mfp = ltl_F(bag_entity("mfp") > 0) f_reach_hspmfp = Formula_rsLTL.f_F(BagDescription.f_TRUE(), (BagDescription.f_entity("hsp:mfp") > 0) ) diff --git a/rsltl_shortcuts.py b/rsltl_shortcuts.py new file mode 100644 index 0000000..815dce4 --- /dev/null +++ b/rsltl_shortcuts.py @@ -0,0 +1,26 @@ +from logics import * + +##### SHORTCUTS + +def bag_True(): + return BagDescription.f_TRUE() + +def bag_entity(name): + return BagDescription.f_entity(name) + +def ltl_F(a0, ctx_arg=None): + if ctx_arg == None: + ctx_arg = bag_True() + return Formula_rsLTL.f_F(ctx_arg, a0) + +def ltl_G(a0, ctx_arg=None): + if ctx_arg == None: + ctx_arg = bag_True() + return Formula_rsLTL.f_G(ctx_arg, a0) + +def ltl_X(a0, ctx_arg=None): + if ctx_arg == None: + ctx_arg = bag_True() + return Formula_rsLTL.f_X(ctx_arg, a0) + + From 32d31f9fa12184df9f56af9c9e763eb9684f892e Mon Sep 17 00:00:00 2001 From: Artur Meski Date: Sat, 2 Dec 2017 00:13:08 +0000 Subject: [PATCH 35/35] Easier formulae --- logics/bags.py | 18 +++-- logics/rsltl.py | 30 ++++++--- logics/rsltl_encoder.py | 12 +++- ...action_system_with_concentrations_param.py | 4 +- rs_testing.py | 35 ++++++++-- rsltl_shortcuts.py | 67 ++++++++++++++++--- 6 files changed, 131 insertions(+), 35 deletions(-) diff --git a/logics/bags.py b/logics/bags.py index b69d694..fd386e7 100644 --- a/logics/bags.py +++ b/logics/bags.py @@ -19,11 +19,11 @@ class BagDescription(object): 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: @@ -37,14 +37,14 @@ class BagDescription(object): def sanity_check(self): """Sanity checks""" - for operand in (self.left_operand, self.right_operand): if operand: if not( isinstance(operand, BagDescription) or isinstance(operand, int)): raise RuntimeError( - "Unexpected operand type for a bag: " + str(operand)) + "Unexpected operand type for a bag: {:s} (type: {:s})".format( + str(operand), str(type(operand)))) @classmethod def f_entity(cls, entity_name): @@ -54,6 +54,14 @@ class BagDescription(object): def f_TRUE(cls): return cls(BagDesc_oper.true) + @classmethod + def f_And(cls, arg_L, arg_R): + return cls(BagDesc_oper.l_and, L_oper=arg_L, R_oper=arg_R) + + @classmethod + def f_Not(cls, arg): + return cls(BagDesc_oper.l_not, L_oper=arg) + def __lt__(self, other): return BagDescription(BagDesc_oper.lt, L_oper=self, R_oper=other) diff --git a/logics/rsltl.py b/logics/rsltl.py index 83d61d1..91fdcc1 100644 --- a/logics/rsltl.py +++ b/logics/rsltl.py @@ -18,6 +18,7 @@ class Formula_rsLTL(object): 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)) @@ -26,34 +27,37 @@ class Formula_rsLTL(object): 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( + return "(" + repr( self.left_operand) + " U[" + repr( self.sub_operand) + "] " + repr( - self.right_operand) + " )" + self.right_operand) + ")" if self.f_type == rsLTL_form_type.t_release: - return "( " + repr( + return "(" + repr( self.left_operand) + " R[" + repr( self.sub_operand) + "] " + repr( - self.right_operand) + " )" + self.right_operand) + ")" @property def is_bag(self): @@ -63,6 +67,10 @@ class Formula_rsLTL(object): def f_bag(cls, bag_descr): return cls(rsLTL_form_type.bag, bag=bag_descr) + @classmethod + def f_Not(cls, arg): + return cls(rsLTL_form_type.l_not, L_oper=arg) + @classmethod def f_And(cls, arg_L, arg_R): return cls(rsLTL_form_type.l_and, L_oper=arg_L, R_oper=arg_R) diff --git a/logics/rsltl_encoder.py b/logics/rsltl_encoder.py index c1f61e7..222535c 100644 --- a/logics/rsltl_encoder.py +++ b/logics/rsltl_encoder.py @@ -120,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): @@ -167,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/reaction_system_with_concentrations_param.py b/rs/reaction_system_with_concentrations_param.py index 38dcf11..4ed6726 100644 --- a/rs/reaction_system_with_concentrations_param.py +++ b/rs/reaction_system_with_concentrations_param.py @@ -23,7 +23,6 @@ class ReactionSystemWithConcentrationsParam(ReactionSystem): self.reactions = [] self.parameters = dict() - # self.parametric_reactions = [] self.meta_reactions = dict() self.permanent_entities = dict() self.background_set = [] @@ -91,8 +90,7 @@ class ReactionSystemWithConcentrationsParam(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_testing.py b/rs_testing.py index 05586bc..5d0e34b 100644 --- a/rs_testing.py +++ b/rs_testing.py @@ -41,10 +41,21 @@ def gene_expression(cmd_args): 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)]) - r.add_reaction([("x",1),("xp",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)]) @@ -66,6 +77,22 @@ def gene_expression(cmd_args): 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(): @@ -89,7 +116,7 @@ def trivial_param(): rc.show() smt_rsc = SmtCheckerRSCParam(rc) - f1 = ltl_F(bag_entity("final") >= 1) + f1 = ltl_F(True, bag_entity("final") >= 1) # f1 = Formula_rsLTL.f_F( # BagDescription.f_TRUE(), @@ -207,7 +234,7 @@ def heat_shock_response_param(cmd_args, print_system=True): # # f_reach_mfp = Formula_rsLTL.f_F(BagDescription.f_TRUE(), (BagDescription.f_entity("mfp") > 0) ) - f_reach_mfp = ltl_F(bag_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) ) diff --git a/rsltl_shortcuts.py b/rsltl_shortcuts.py index 815dce4..1e3abcc 100644 --- a/rsltl_shortcuts.py +++ b/rsltl_shortcuts.py @@ -8,19 +8,68 @@ def bag_True(): def bag_entity(name): return BagDescription.f_entity(name) -def ltl_F(a0, ctx_arg=None): - if ctx_arg == None: - ctx_arg = bag_True() +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(a0, ctx_arg=None): - if ctx_arg == None: - ctx_arg = bag_True() +def ltl_G(ctx_arg, a0): + a0 = get_bag_if_str(a0) return Formula_rsLTL.f_G(ctx_arg, a0) -def ltl_X(a0, ctx_arg=None): - if ctx_arg == None: - ctx_arg = bag_True() +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)