Cleanup, example

rs_testing.py

@@ -26,6 +26,80 @@ def gene_expression(cmd_args):
     """
     r = ReactionSystemWithConcentrationsParam()
     
+    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(("h", 1))
+
+    all_entities = set(r.background_set)
+
+    ## y
+    # 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)])
+    lda1 = r.get_param("lda1")
+    lda2 = r.get_param("lda2")
+    lda3 = r.get_param("lda3")
+
+    r.add_reaction([("y",1)],[("h",1)],lda1)
+    r.add_reaction(lda2,[("h",1)],[("yp",1)])
+    r.add_reaction([("y",1),("yp",1)],[("h",1)],lda3)
+
+    c = ContextAutomatonWithConcentrations(r)
+    c.add_init_state("0")
+    c.add_state("1")
+    
+    # the experiments starts with adding x and y:
+    c.add_transition("0", [("y", 1)], "1")
+    
+    # for all the remaining steps we have empty context sequences
+    c.add_transition("1", [], "1")
+
+    rc = ReactionSystemWithAutomaton(r, c)
+    rc.show()
+    
+    # f_x1 = ltl_G(True, ltl_Implies(
+    #     exact_state(["x"], all_entities),
+    #     ltl_X(bag_And(bag_Not("Y"), bag_Not("Z"), bag_Not("h")), exact_state(["x", "xp"], all_entities))))
+    # f_x2 = ltl_G(True, ltl_Implies(
+    #     exact_state(["x", "xp"], all_entities),
+    #     ltl_X(bag_Not("h"), exact_state("X", all_entities))))
+
+    reach_yp = ltl_F(True, "yp")
+    reach_Y = ltl_F(True, "Y")
+    reach_y = ltl_F(True, "y")
+
+    f_y1 = ltl_G(bag_entity("h") == 0, ltl_Implies(
+        bag_entity("y") > 0,
+        ltl_X(True, bag_And(bag_entity("y") > 0, bag_entity("yp") > 0))))
+
+    f_y2 = ltl_G(bag_entity("h") == 0, ltl_Implies(
+        bag_And(bag_entity("y") > 0, bag_entity("yp") > 0),
+        ltl_F(True, bag_entity("Y") > 0)))
+        
+    # delayed_Y = ltl_X(True, ltl_And(bag_entity("Y") == 0, ltl_X(True, bag_entity("Y") == 0)))
+    delayed_Y = ltl_X(True, bag_entity("Y") == 0)
+
+    obs_1 = ltl_And(f_y1, reach_y, delayed_Y)
+    obs_2 = ltl_And(f_y2, reach_y, reach_Y, delayed_Y)
+
+    # f_x2 = ltl_G(True, ltl_Implies(
+    #     exact_state(["x", "xp"], all_entities),
+    #     ltl_X(bag_Not("h"), exact_state("X", all_entities))))
+    #
+        
+    smt_rsc = SmtCheckerRSCParam(rc, optimise=cmd_args.optimise)
+    #
+    smt_rsc.check_rsltl(formulae_list=[obs_1, obs_2], max_level=4, cont_if_sat=True)
+
+    # smt_rsc.check_rsltl(formula=f_x1, print_witness=True)
+
+def gene_expression_full(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))
@@ -45,23 +119,24 @@ def gene_expression(cmd_args):
     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)])
+    ## y
-    param_p1 = r.get_param("P1")
+    # r.add_reaction([("y",1)],[("h",1)],[("y",1)])
-    r.add_reaction(param_p1,[("h",1)],[("xp",1)])
+    # r.add_reaction([("y",1)],[("Q",1)],[("yp",1)])
-    
+    # r.add_reaction([("y",1),("yp",1)],[("h",1)],[("Y",1)])
-    # r.add_reaction([("x",1),("xp",1)],[("h",1)],[("X",1)])
+    lda1 = r.get_param("lda1")
-    param_p2_r = r.get_param("P2_r")
+    lda2 = r.get_param("lda2")
-    param_p2_i = r.get_param("P2_i")
+    lda3 = r.get_param("lda3")
-    param_p2_p = r.get_param("P2_p")
+    r.add_reaction([("y",1)],[("h",1)],lda1)
-    r.add_reaction(param_p2_r,param_p2_i,param_p2_p)
+    r.add_reaction(lda2,[("Q",1)],[("yp",1)])
+    r.add_reaction([("y",1),("yp",1)],[("h",1)],lda3)
 
-    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)
@@ -77,19 +152,36 @@ def gene_expression(cmd_args):
     rc = ReactionSystemWithAutomaton(r, c)
     rc.show()
     
-    f_x1 = ltl_G(True, ltl_Implies(
+    # f_x1 = ltl_G(True, ltl_Implies(
-        exact_state(["x"], all_entities),
+    #     exact_state(["x"], all_entities),
-        ltl_X(bag_And(bag_Not("Y"), bag_Not("Z"), bag_Not("h")), exact_state(["x", "xp"], 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(
+    # f_x2 = ltl_G(True, ltl_Implies(
-        exact_state(["x", "xp"], all_entities),
+    #     exact_state(["x", "xp"], all_entities),
-        ltl_X(bag_Not("h"), exact_state("X", all_entities))))
+    #     ltl_X(bag_Not("h"), exact_state("X", all_entities))))
-    
+
-    reach_xp = ltl_F(True, "xp")
+    reach_yp = ltl_F(True, "yp")
-    reach_X = ltl_F(True, "X")
+    reach_Y = ltl_F(True, "Y")
+    reach_y = ltl_F(True, "y")
+
+    f_y1 = ltl_G(bag_entity("h") == 0, ltl_Implies(
+        bag_And(bag_entity("Q") == 0, bag_entity("y") > 0),
+        ltl_X(True, bag_And(bag_entity("y") > 0, bag_entity("yp") > 0))))
+
+    f_y2 = ltl_G(bag_entity("h") == 0, ltl_Implies(
+        bag_And(bag_entity("y") > 0, bag_entity("yp") > 0),
+        ltl_F(True, bag_entity("Y") > 0)))
+
+    obs_1 = ltl_And(f_y1, reach_y)
+    obs_2 = ltl_And(f_y2, reach_y, reach_Y)
+
+    # f_x2 = ltl_G(True, ltl_Implies(
+    #     exact_state(["x", "xp"], all_entities),
+    #     ltl_X(bag_Not("h"), exact_state("X", all_entities))))
+    #
         
     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(formulae_list=[obs_1, obs_2], max_level=3, cont_if_sat=True)
 
     # smt_rsc.check_rsltl(formula=f_x1, print_witness=True)
     

rsltl_shortcuts.py

@@ -27,6 +27,10 @@ def bag_And(*args):
     return last
     
 def exact_state(contained_entities, all_entities):
+    """
+    Assumes 0 concentration level for all the
+    not listed entities but present in all_entities
+    """
     expr = []
     for ent_str in all_entities:
         ent = bag_entity(ent_str)

smt/smt_checker_rsc_param.py

@@ -161,8 +161,8 @@ class SmtCheckerRSCParam(object):
         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[path_idx]) == level
-        assert len(self.path_v_improd_for_entities[path_idx]) == level
+        # assert len(self.path_v_improd_for_entities[path_idx]) == level
 
         if level < 1:
             #
@@ -172,45 +172,44 @@ class SmtCheckerRSCParam(object):
             #
             self.path_v_improd[path_idx].append(None)
             self.path_v_improd_for_entities[path_idx].append(None)
-        else:
 
-            reactions_dict = dict()
+        reactions_dict = dict()
-            number_of_reactions = len(self.rs.reactions)
+        number_of_reactions = len(self.rs.reactions)
 
-            all_entities_dict = dict()
+        all_entities_dict = dict()
 
-            for reaction in self.rs.reactions:
+        for reaction in self.rs.reactions:
-                *_, products = reaction
+            *_, products = reaction
 
-                reaction_id = self.rs.reactions.index(reaction)
+            reaction_id = self.rs.reactions.index(reaction)
 
-                entities_dict = dict()
+            entities_dict = dict()
-                
-                if is_param(products):
-
-                    for entity in self.rs.set_of_bgset_ids:
-                        entity_name = self.rs.get_entity_name(entity)
-                        new_var = Int("p{:d}L{:d}_ImProd_r{:d}_{:s}".format(path_idx, level, reaction_id, entity_name))
-                        entities_dict[entity] = new_var
-
-                        all_entities_dict.setdefault(entity, [])
-                        all_entities_dict[entity].append(new_var)
-                                            
-                else:
-                
-                    for entity, conc in products:
-                        entity_name = self.rs.get_entity_name(entity)
-                        new_var = Int("p{:d}L{:d}_ImProd_r{:d}_{:s}".format(path_idx, level, reaction_id, entity_name))
-                        entities_dict[entity] = new_var
-
-                        all_entities_dict.setdefault(entity, [])
-                        all_entities_dict[entity].append(new_var)
-
-                reactions_dict[reaction_id] = entities_dict
             
-            self.path_v_improd[path_idx].append(reactions_dict)
+            if is_param(products):
-            self.path_v_improd_for_entities[path_idx].append(all_entities_dict)
+
+                for entity in self.rs.set_of_bgset_ids:
+                    entity_name = self.rs.get_entity_name(entity)
+                    new_var = Int("p{:d}L{:d}_ImProd_r{:d}_{:s}".format(path_idx, level, reaction_id, entity_name))
+                    entities_dict[entity] = new_var
+
+                    all_entities_dict.setdefault(entity, [])
+                    all_entities_dict[entity].append(new_var)
+                                        
+            else:
+            
+                for entity, conc in products:
+                    entity_name = self.rs.get_entity_name(entity)
+                    new_var = Int("p{:d}L{:d}_ImProd_r{:d}_{:s}".format(path_idx, level, reaction_id, entity_name))
+                    entities_dict[entity] = new_var
+
+                    all_entities_dict.setdefault(entity, [])
+                    all_entities_dict[entity].append(new_var)
+
+            reactions_dict[reaction_id] = entities_dict
         
+        self.path_v_improd[path_idx].append(reactions_dict)
+        self.path_v_improd_for_entities[path_idx].append(all_entities_dict)
+    
     def prepare_param_variables(self):
         """
         Prepares variables for parameters
@@ -270,8 +269,9 @@ class SmtCheckerRSCParam(object):
 
         for param_vars in self.v_param.values():
             for pvar in param_vars:
-                self.solver.add_soft(pvar < 1)
+                # self.solver.add_soft(pvar == 0)
-                                        
+                self.solver.minimize(pvar)
+                                       
     def enc_concentration_levels_assertion(self, level, path_idx):
         """
         Encodes assertions that (some) variables need to be >=0
@@ -599,10 +599,48 @@ class SmtCheckerRSCParam(object):
                             end="")
                 print(" }")
         
-        # Parameters    
+        print()
         
-        print("\n\n  Parameters:\n")
+    def get_enc_formulae(self, encoder, formulae_list):
-        for param_name in self.rs.parameters.keys():
+        enc_form = []
+        for formula in formulae_list:
+            
+            path_idx = formulae_list.index(formula)
+            
+            print_info("Generating the encoding for {:s} ({:d} of {:d})".format(
+                str(formula), path_idx+1, len(formulae_list)))
+
+            encoder.load_variables(
+                var_rs=self.path_v[path_idx],
+                var_ctx=self.path_v_ctx[path_idx],
+                var_loop_pos=self.path_loop_position[path_idx])
+                
+            enc_form.append(encoder.get_encoding(formula, self.current_level))
+            ncalls = encoder.get_ncalls()
+
+            print_info("Cache hits: {:d}, encode calls: {:d} (approx: {:d})".format(
+                    encoder.get_cache_hits(), ncalls[0], ncalls[1]))
+
+            encoder.flush_cache()
+            
+        return enc_form
+        
+    def print_witness(self, formulae_list):
+    
+        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)
+        
+        self.print_parameter_valuations()   
+    
+    def print_parameter_valuations(self):
+
+        m = self.solver.model()
+
+        print("\n  Parameters:\n")
+        for param_name in sorted(self.rs.parameters.keys()):
             print("{: >6}: ".format(param_name), end="")
             print("{", end="")
             
@@ -618,9 +656,28 @@ class SmtCheckerRSCParam(object):
                         " " + str(self.rs.get_entity_name(entity)) + "=" + str(var_rep),
                         end="")
             print(" }")
+    
+        print()
         
-        print("\n")
+    def enc_concentration_levels_assertions_for_paths(self, level, num_of_paths):
-            
+    
+        additional_assertions = []
+        for path_idx in range(num_of_paths):
+            additional_assertions.append(self.enc_concentration_levels_assertion(level, path_idx))
+        additional_assertions.append(self.enc_param_concentration_levels_assertion())
+        
+        return additional_assertions
+    
+    def enc_transition_relation_for_paths(self, level, num_of_paths):
+        enc_trans = []
+        for path_idx in range(num_of_paths):
+            enc_trans.append(self.enc_transition_relation(level, path_idx))
+        return enc_trans
+        
+    def print_level(self):
+        print(
+            "{:->70}".format("[ level=" + str(self.current_level) + " done ]"))
+    
     def check_rsltl(
             self, formulae_list, 
             print_witness=True, 
@@ -653,10 +710,12 @@ class SmtCheckerRSCParam(object):
         for form in formulae_list:
             print_info(" "*4 + str(form))
 
-        if print_time:
+        print_info("INITIALISING...")
-            start = resource.getrusage(resource.RUSAGE_SELF).ru_utime
 
-        self.prepare_all_variables(num_of_paths)        
+        if print_time:
+            start_time = resource.getrusage(resource.RUSAGE_SELF).ru_utime
+
+        self.prepare_all_variables(num_of_paths)
         
         self.load_varset_for_path(0)
         
@@ -668,26 +727,21 @@ class SmtCheckerRSCParam(object):
         
         self.current_level = 0
 
-        self.prepare_all_variables(num_of_paths)
+        # self.prepare_all_variables(num_of_paths)
 
         # assertions for all the paths and parameters
-        additional_assertions = []
+        self.solver_add(self.enc_concentration_levels_assertions_for_paths(0, num_of_paths))
-        for path_idx in range(num_of_paths):
+        self.solver_add(self.enc_param_concentration_levels_assertion())
-            additional_assertions.append(self.enc_concentration_levels_assertion(0, path_idx))
+
-        additional_assertions.append(self.enc_param_concentration_levels_assertion())
+        encoder = rsLTL_Encoder(self)
-        self.solver_add(additional_assertions)
 
         if self.optimise:
             self.assert_param_optimisation()
 
-        encoder = rsLTL_Encoder(self)
+        print_info("STARTING TO ITERATE...")
-
-        print_info("Iterating...")
         
         while True:
             
-
-            
             print(
                 "\n{:-^70}".format("[ Working at level=" + str(self.current_level) + " ]"))
             # stdout.flush()
@@ -695,36 +749,18 @@ class SmtCheckerRSCParam(object):
             # reachability test:
             self.solver.push()
 
-            enc_form = []
+            # Encoding all the formulae
-            for formula in formulae_list:
+            enc_form = self.get_enc_formulae(encoder, formulae_list)
-                
-                path_idx = formulae_list.index(formula)
-                
-                print_info("Generating the encoding for {:s} ({:d} of {:d})".format(
-                    str(formula), path_idx+1, len(formulae_list)))
-
-                encoder.load_variables(
-                    var_rs=self.path_v[path_idx],
-                    var_ctx=self.path_v_ctx[path_idx],
-                    var_loop_pos=self.path_loop_position[path_idx])
-                    
-                enc_form.append(encoder.get_encoding(formula, self.current_level))
-                ncalls = encoder.get_ncalls()
-
-                print_info("Cache hits: {:d}, encode calls: {:d} (approx: {:d})".format(
-                        encoder.get_cache_hits(), ncalls[0], ncalls[1]))
-
-                encoder.flush_cache()
-                
             print_info("Adding the formulae to the solver...")
-            
-            # print (enc_form)
-            
             self.solver_add(enc_form)
 
-            print_info("Adding the loops encoding...")
+            # Loops encoding
+            print_info("Adding the encoding for the loops...")
             self.solver_add(self.get_loop_encodings())
 
+            # if self.optimise:
+            #     self.assert_param_optimisation()
+
             print_info("Testing satisfiability...")
             result = self.solver.check()
             if result == sat:
@@ -732,55 +768,49 @@ class SmtCheckerRSCParam(object):
                         "SAT at level={:d}".format(self.current_level)))
                 # print(self.solver.model())
                 if print_witness:
-                    for formula in formulae_list:
+                    self.print_witness(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")
+                print_info("UNSAT")
                 self.solver.pop()
 
             self.prepare_all_variables(num_of_paths)
             
             # assertions for all the paths
-            additional_assertions = []
+            self.solver_add(self.enc_concentration_levels_assertions_for_paths(self.current_level + 1, num_of_paths))
-            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_for_paths(self.current_level, num_of_paths))
-                self.solver_add(self.enc_transition_relation(self.current_level, path_idx))
+
+            self.print_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))
+                print_info("Stopping at level=" + str(max_level))
                 break
 
         if print_time:
-            # stop = time()
+            self.print_time(start_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(
+            self.print_mem()
-                "[i] {: >60}".format(
+
-                    " Memory: " +
+    def print_time(self, start):
-                    repr(
+        stop = resource.getrusage(resource.RUSAGE_SELF).ru_utime
-                        resource.getrusage(resource.RUSAGE_SELF).ru_maxrss /
+        self.verification_time = stop - start
-                        (1024 * 1024)) + " MB"))
+        print()
-        
+        print_info("{: >60}".format(
+                " Time: " + repr(self.verification_time) + " s"))
+    
+    def print_mem(self):
+        print_info(
+            "{: >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"""