changes for emsoft 2016 camera ready (gearbox runs)

hylaa/aggdag.py

@@ -98,7 +98,7 @@ def add_invariant_op(self, step, i_index, is_stronger):
         self.cur_node.op_list.append(op)
 
     def make_op_transition(self, t, t_lpi, state, parent_node):
-        'make an OpTransition object'
+        'make an OpTransition object, can return null if lp solving fails'
 
         step_in_mode = state.cur_step_in_mode
         steps_since_start = state.cur_steps_since_start
@@ -108,19 +108,22 @@ def make_op_transition(self, t, t_lpi, state, parent_node):
 
         op = OpTransition(step_in_mode, parent_node, None, t, None)
 
-        self.settings.aggstrat.pretransition(t, t_lpi, op)
+        succeeded = self.settings.aggstrat.pretransition(t, t_lpi, op)
 
-        lputil.add_reset_variables(t_lpi, t.to_mode.mode_id, t.transition_index, \
-            reset_csr=t.reset_csr, minkowski_csr=t.reset_minkowski_csr, \
-            minkowski_constraints_csr=t.reset_minkowski_constraints_csr, \
-            minkowski_constraints_rhs=t.reset_minkowski_constraints_rhs, successor_has_inputs=successor_has_inputs)
+        if not succeeded:
+            op = None
+        else:
+            lputil.add_reset_variables(t_lpi, t.to_mode.mode_id, t.transition_index, \
+                reset_csr=t.reset_csr, minkowski_csr=t.reset_minkowski_csr, \
+                minkowski_constraints_csr=t.reset_minkowski_constraints_csr, \
+                minkowski_constraints_rhs=t.reset_minkowski_constraints_rhs, successor_has_inputs=successor_has_inputs)
 
-        if not t_lpi.is_feasible():
-            raise RuntimeError("cur_state became infeasible after reset was applied")
+            if not t_lpi.is_feasible():
+                raise RuntimeError("cur_state became infeasible after reset was applied")
 
-        op_list = [op]
-        state = StateSet(t_lpi, t.to_mode, steps_since_start, op_list, is_concrete)
-        op.poststate = state
+            op_list = [op]
+            state = StateSet(t_lpi, t.to_mode, steps_since_start, op_list, is_concrete)
+            op.poststate = state
 
         return op
 
@@ -147,12 +150,17 @@ def _get_node_leaf_ops(self, node):
 
         rv = []
 
+        child_nodes = []
+
         for op in node.op_list:
             if isinstance(op, OpTransition):
                 if op.child_node is None:
                     rv.append(op)
-                else:
-                    rv += self._get_node_leaf_ops(op.child_node)
+                elif not op.child_node in child_nodes: # don't insert the same node twice
+                    child_nodes.append(op.child_node)
+
+        for child_node in child_nodes:
+            rv += self._get_node_leaf_ops(child_node)
 
         return rv
 
@@ -424,6 +432,7 @@ def replay_op(self, op_list, i):
             if True or not skip:
                 self.aggdag.core.print_verbose(
                     f"doing invariant intersection in replay at step {cur_state.cur_step_in_mode}")
+
                 is_feasible = self.replay_op_intersect_invariant(cur_state, op)
 
                 if not is_feasible:
@@ -453,17 +462,21 @@ def replay_op_transition(self, state, op):
                 self.aggdag.core.error_reached(state, t, t_lpi)
 
             new_op = self.aggdag.make_op_transition(t, t_lpi, state, self)
-            self.op_list.append(new_op)
 
-            if op.child_node is None:
-                self.aggdag.waiting_list.append(new_op)
-                print_verbose("Replay Transition {} when deaggreaged to steps {}".format( \
-                              t, state.cur_steps_since_start))
+            if not new_op: # make_op_transition can fail due to numerical issues
+                print_verbose("Replay Transition {} became infeasible after changing opt direction, skipping")
             else:
-                self.aggdag.deagg_man.update_transition_successors(op, new_op)
-
-                print_verbose("Replay Transition refined transition {} when deaggregated at steps {}".format( \
-                              t, state.cur_steps_since_start))
+                self.op_list.append(new_op)
+
+                if op.child_node is None:
+                    self.aggdag.waiting_list.append(new_op)
+                    print_verbose("Replay Transition {} when deaggreaged to steps {}".format( \
+                                  t, state.cur_steps_since_start))
+                else:
+                    self.aggdag.deagg_man.update_transition_successors(op, new_op)
+
+                    print_verbose("Replay Transition refined transition {} when deaggregated at steps {}".format( \
+                                  t, state.cur_steps_since_start))
         else:
             print_verbose(f"Replay skipped transition {t} when deaggregated to steps {state.cur_steps_since_start}")
 
@@ -478,8 +491,14 @@ def replay_op_intersect_invariant(self, state, op):
         state.step(step)
 
         lc = state.mode.inv_list[invariant_index]
+
+        has_intersection = lputil.check_intersection(state.lpi, lc.negate())
 
-        if lputil.check_intersection(state.lpi, lc.negate()):
+        if has_intersection is None:
+            rv = False # not feasible
+        elif not has_intersection:
+            rv = True # no intersection and still feasible
+        else: 
             old_row = state.invariant_constraint_rows[invariant_index]
             vec = lc.csr.toarray()[0]
             rhs = lc.rhs
@@ -500,7 +519,9 @@ def replay_op_intersect_invariant(self, state, op):
             op = OpInvIntersect(step, self, invariant_index, is_stronger)
             self.op_list.append(op)
 
-        return state.lpi.is_feasible()
+            rv = state.lpi.is_feasible()
+
+        return rv
 
     def get_cur_state(self):
         'get the current state for this node (None if it has left the invariant)'

hylaa/aggstrat.py

@@ -177,9 +177,16 @@ def get_simulation_pop_mode(self, sim_waiting_list):
     def pretransition(self, t, t_lpi, op_transition):
         'event function, called when taking a transition before the reset is applied'
 
+        rv = True
+
         if self.agg_type == Aggregated.AGG_ARNOLDI_BOX:
             op_transition.premode_center = lputil.get_box_center(t_lpi)
 
+            if op_transition.premode_center is None:
+                rv = False
+
+        return rv
+
     def get_agg_type(self, op_list):
         '''
         Gets the type of aggregation to be performed for the passed in objects. 

hylaa/core.py

@@ -22,7 +22,7 @@
 class Core(Freezable):
     'main computation object. initialize and call run()'
 
-    def __init__(self, ha, hylaa_settings, seed=0):
+    def __init__(self, ha, hylaa_settings):
         assert isinstance(hylaa_settings, HylaaSettings)
         assert isinstance(ha, HybridAutomaton)
 
@@ -52,7 +52,7 @@ def __init__(self, ha, hylaa_settings, seed=0):
         self.sim_took_transition = None  # list of booleans for each sim
 
         # make random number generation (for example, to find orthogonal directions) deterministic
-        np.random.seed(seed=seed)
+        np.random.seed(hylaa_settings.random_seed)
 
         LpInstance.print_normal = self.print_normal
         LpInstance.print_verbose = self.print_verbose
@@ -447,6 +447,9 @@ def simulate(self, init_mode, box, num_sims):
         # initialize time elapse in each mode of the hybrid automaton
         ha = init_mode.ha
 
+        # initialize result object (many fields unused during simulation)
+        self.result = HylaaResult()
+
         self.setup_ha(ha)
 
         # each simulation is a tuple (mode, pt, num_steps)
@@ -471,6 +474,8 @@ def simulate(self, init_mode, box, num_sims):
 
         self.plotman.compute_and_animate()
 
+        return self.result
+
     def sim_pop_waiting_list(self):
         'pop a state off the simulation waiting list'
 
@@ -501,13 +506,14 @@ def do_step_sim(self):
         'do a simulation step'
 
         if not self.sim_states:
-            # pop minimum time state
-            self.sim_pop_waiting_list()
-
-            self.sim_should_try_guards = self.settings.process_urgent_guards
-            self.sim_took_transition = [False] * len(self.sim_states)
-            self.plotman.pause()
-            self.print_verbose("Pausing due to sim pop()")
+            if self.sim_waiting_list: # if simulation is not done
+                # pop minimum time state
+                self.sim_pop_waiting_list()
+
+                self.sim_should_try_guards = self.settings.process_urgent_guards
+                self.sim_took_transition = [False] * len(self.sim_states)
+                self.plotman.pause()
+                self.print_verbose("Pausing due to sim pop()")
         else:
             # simulate one step for all states in sim_states
             finished = True

hylaa/deaggregation.py

@@ -8,6 +8,7 @@
 from collections import deque, namedtuple
 
 from hylaa.util import Freezable
+from hylaa.timerutil import Timers
 
 # Operation types
 OpInvIntersect = namedtuple('OpInvIntersect', ['step', 'node', 'i_index', 'is_stronger'])
@@ -193,8 +194,12 @@ def update_transition_successors(self, old_op, new_op):
     def begin_replay(self, node):
         'begin a deaggregation replay with the passed-in node'
 
+        Timers.tic('begin deagg replay')
+
         # remove all states in the waiting list that come from this node
         self.aggdag.remove_node_decendants_from_waiting_list(node)
 
         # start to populate waiting_nodes
         self.waiting_nodes.append((node, node.split()))
+
+        Timers.toc('begin deagg replay')

hylaa/hybrid_automaton.py

@@ -427,12 +427,15 @@ def get_guard_intersection(self, lpi):
 
             lpi.set_minimize_direction(row, is_csr=True)
 
-            #print('.1 hybrid_automaton t={}, is_feasible before = {}'.format(i, lpi.is_feasible()))
-            #print('.2 hybrid_automaton t={}, is_feasible before = {}'.format(i, lpi.is_feasible()))
-
             columns = [lpi.cur_vars_offset + i for i in row.indices]
 
-            result = lpi.minimize(columns=columns)
+            result = lpi.minimize(columns=columns, fail_on_unsat=False)
+
+            # sometimes, changing the objective function makes lp infeasible (due to numerical precision issues)
+            # this happens on gearbox with small time steps. in this case, just return no transition is possible
+            if result is None:
+                all_sat = False
+                break
 
             dot_res = np.dot(result, row.data)
 

hylaa/lpinstance.py

@@ -115,7 +115,7 @@ def _create_bm_indices(self):
     def _column_names_str(self, cur_var_print):
         'get the line in __str__ for the column names'
 
-        rv = "   "
+        rv = "       "
 
         for col, name in enumerate(self.names):
             name = self.names[col]
@@ -140,7 +140,7 @@ def _opt_dir_str(self, zero_print):
 
         lp = self.lp
         cols = self.get_num_cols()
-        rv = "min"
+        rv = "min    "
 
         for col in range(1, cols + 1):
             val = glpk.glp_get_obj_coef(lp, col)
@@ -537,7 +537,7 @@ def set_minimize_direction(self, direction_vec, is_csr=False, offset=None):
 
         Timers.toc("set_minimize_direction")
 
-    def minimize(self, direction_vec=None, columns=None, fail_on_unsat=True):
+    def minimize(self, direction_vec=None, columns=None, fail_on_unsat=True, print_on=False):
         '''minimize the lp, returning a list of assigments to each of the variables
 
         if direction_vec is not None, this will first assign the optimization direction (note: relative to cur_vars)
@@ -557,7 +557,7 @@ def minimize(self, direction_vec=None, columns=None, fail_on_unsat=True):
         params = glpk.glp_smcp()
         glpk.glp_init_smcp(params)
         params.meth = glpk.GLP_DUALP # use dual simplex since we're reoptimizing often
-        params.msg_lev = glpk.GLP_MSG_OFF
+        params.msg_lev = glpk.GLP_MSG_ALL if print_on else glpk.GLP_MSG_OFF
         params.tm_lim = 1000 # 1000 ms time limit
 
         Timers.tic('glp_simplex')
@@ -586,10 +586,10 @@ def minimize(self, direction_vec=None, columns=None, fail_on_unsat=True):
 
         if rv is None and fail_on_unsat:
             LpInstance.print_normal("Note: minimize failed with fail_on_unsat was true, resetting and retrying...")
-                        
+
             glpk.glp_cpx_basis(self.lp) # resets the initial basis
 
-            rv = self.minimize(direction_vec, columns, False)
+            rv = self.minimize(direction_vec, columns, False, print_on=True)
 
             if rv is not None:
                 LpInstance.print_verbose("Note: LP was infeasible, but then feasible after resetting statuses")

hylaa/lputil.py

@@ -304,6 +304,10 @@ def check_intersection(lpi, lc, tol=1e-13):
     '''check if there is an intersection between the LP constriants and the LinearConstraint object lc
 
     This solves an LP optimizing in the given direction... without adding the constraint to the LP
+
+    This returns True/False if an intersection is possible
+
+    it also can return None if the lp is infeasible
     '''
 
     Timers.tic("check_intersection")
@@ -313,13 +317,19 @@ def check_intersection(lpi, lc, tol=1e-13):
     columns = lc.csr.indices[0:lc.csr.indptr[1]]
     lp_columns = [lpi.cur_vars_offset + c for c in columns]
 
-    lp_res = lpi.minimize(columns=lp_columns)
+    lp_res = lpi.minimize(columns=lp_columns, fail_on_unsat=False)
 
-    dot_res = np.dot(lc.csr.data, lp_res)
+    if lp_res is None:
+        # sometimes, changing optimization direction makes lp infeasible (up to numerical accuracy)
+        # this happens in gearbox with small time steps. In this case, return no intersection
+        rv = None
+    else:
+        dot_res = np.dot(lc.csr.data, lp_res)
+        rv = dot_res + tol <= lc.rhs
 
     Timers.toc("check_intersection")
 
-    return dot_res + tol <= lc.rhs
+    return rv
 
 def add_init_constraint(lpi, vec, rhs, basis_matrix=None, input_effects_list=None, row_index=None):
     '''
@@ -831,7 +841,10 @@ def add_curtime_constraints(lpi, csr, rhs_vec):
     lpi.set_constraints_csr(csr, offset=(prerows, lpi.cur_vars_offset))
 
 def get_box_center(lpi):
-    '''get the center of the box overapproximation of the passed-in lpi'''
+    '''get the center of the box overapproximation of the passed-in lpi
+
+    may return None if lp solving fails (numerical issues)
+    '''
 
     Timers.tic('get_box_center')
 
@@ -843,8 +856,15 @@ def get_box_center(lpi):
         min_dir = [1 if i == dim else 0 for i in range(dims)]
         max_dir = [-1 if i == dim else 0 for i in range(dims)]
 
-        min_val = lpi.minimize(direction_vec=min_dir, columns=[col])[0]
-        max_val = lpi.minimize(direction_vec=max_dir, columns=[col])[0]
+        min_val = lpi.minimize(direction_vec=min_dir, columns=[col], fail_on_unsat=False)
+        max_val = lpi.minimize(direction_vec=max_dir, columns=[col], fail_on_unsat=False)
+
+        if min_val is None or max_val is None:
+            pt = None
+            break
+
+        min_val = min_val[0]
+        max_val = max_val[0]
 
         pt.append((min_val + max_val) / 2.0)