From 9297a7caf7f1e33a655d240096bb501529012fcf Mon Sep 17 00:00:00 2001 From: Artur Meski Date: Sun, 17 Sep 2017 12:24:07 +0100 Subject: [PATCH] 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)