diff --git a/README.md b/README.md index 855a441..e0663d8 100644 --- a/README.md +++ b/README.md @@ -23,6 +23,21 @@ From there, it's a priority first search, prioritizing first: It creates a partial state graph and seeks $K$ satisfying states. +## Installation + +It is recommended to use a virtual environment to locally install the dependencies (since StormPy is a bit of a beast). In the project directory: + +```sh +# Create a virtual environment +python3 -m venv .venv +# Activate the virtual environment +source .venv/bin/activate +# Install the dependencies +pip -r install dependencies.txt +``` + +When done, the `deactivate` command can be used. + ## Usage ```bash diff --git a/all_test.sh b/all_test.sh new file mode 100755 index 0000000..13b10cf --- /dev/null +++ b/all_test.sh @@ -0,0 +1,31 @@ +#!/bin/sh + +# This test script is really hacky. Please don't judge me. + +rm -rf results/isr +rm -rf results/isr_piped +# rm -rf results/single_order +# rm -rf results/single_order_piped + +mkdir -p results/isr +mkdir -p results/isr_piped +# mkdir -p results/single_order +# mkdir -p results/single_order_piped + +# Test ISR +python3 test.py +mv results/*.csv results/isr + +# Test ISR with Piped space scaling +python3 test.py --piped +mv results/*.csv results/isr_piped + +echo "WARNING! Not testing SOP" + +# Test SOP +# python3 test.py --sop +# mv results/*.csv results/single_order +# +# # Test SOP with Piped space scaling +# python3 test.py --sop --piped +# mv results/*.csv results/single_order_piped diff --git a/dependencies.txt b/dependencies.txt new file mode 100644 index 0000000..1f5cebe --- /dev/null +++ b/dependencies.txt @@ -0,0 +1,2 @@ +numpy +stormpy diff --git a/src/crn.py b/src/crn.py index cb916e3..57b090f 100644 --- a/src/crn.py +++ b/src/crn.py @@ -23,7 +23,7 @@ class Transition: # vector #: np.vector # enabled #: lambda # rate_finder #: lambda - def __init__(self, vector, enabled, rate_finder, name=None, rate_constant=None): + def __init__(self, vector, enabled, rate_finder, name=None, rate_constant=None, dim_bounds=None): self.vector = np.array(vector) self.vec_as_mat = np.matrix(vector).T self.enabled_lambda = enabled @@ -31,13 +31,42 @@ def __init__(self, vector, enabled, rate_finder, name=None, rate_constant=None): self.can_eliminate = False self.name = name self.rate_constant = rate_constant + self.in_s0 = False + self.dim_bounds = dim_bounds def enabled(self, state): - return self.enabled_lambda(state) + if self.dim_bounds is None: + return self.enabled_lambda(state) + else: + # This checks to see if we impose an upper limit on the variable at index i and then if so, + # if we are below that upper limit. + idx_passes = lambda i : self.dim_bounds[i] == -1 or self.vector[i] + state[i] < self.dim_bounds[i] + return np.all([idx_passes(i) for i in range(len(state))]) and self.enabled_lambda(state) + +class SortableTransition: + def __init__(self, transition : Transition, index : int): + self.transition = transition + self.rate_constant = transition.rate_constant + self.index = index def __str__(self): return self.name + def __gt__(self, other): + return self.rate_constant > other.rate_constant + + def __le__(self, other): + return self.rate_constant <= other.rate_constant + + def __eq__(self, other): + return self.rate_constant == other.rate_constant + + def __lt__(self, other): + return self.rate_constant < other.rate_constant + + def __ge__(self, other): + return self.rate_constant >= other.rate_constant + class Crn: def __init__(self, transitions=None, boundary=None, init_state=None, all_trans_always_enabled=False, all_rates_const=False): self.boundary = boundary diff --git a/src/cycle.py b/src/cycle.py new file mode 100644 index 0000000..5477bf9 --- /dev/null +++ b/src/cycle.py @@ -0,0 +1,159 @@ +import numpy as np +from scipy.linalg import null_space + +from crn import * +from subspace import * + +# TODO: are we sure we need both? +from fractions import Fraction + +class Cycle: + def __init__(self, ordered_reactions : list): + ''' + Creates an object + ''' + self.__ordered_reactions = ordered_reactions + self.__check_cycle_valid() + self.__in_s0 = [r.in_s0 for r in ordered_reactions] + if not np.any(self.__in_s0): + raise Exception("Cycle must leave S0 (else it may be useless)!") + self.is_orthocycle = np.all(self.__in_s0) + + def __check_cycle_valid(self) -> bool: + sum = self.__ordered_reactions[0].vec_as_mat + for r in self.__ordered_reactions[1::]: + sum += r.vec_as_mat + assert(np.all(np.is_close(sum, 0))) + + def apply_cycle(self, state : np.matrix) -> list: + s = state + new_states = [] + for r in self.__ordered_reactions: + s += r.vec_as_mat + new_states.append(s) + return new_states + +def get_commutable_transitions(crn : Crn, s0 : Subspace, ss : Subspace) -> list: + transitions = [] + for t in crn.transitions: + # TODO: need offset? + if np.all(np.is_close(s0.P * t.vec_as_mat, 0)).all() and np.all(np.isclose(ss.P * t.vec_as_mat, t.vec_as_mat)): + transitions.append(t) + return t + +primes = [2] + +def is_prime(i : int) -> bool: + for j in range(2, i): + if i % j == 0: + return False + return True + +def init_primes_to(n : int): + last_prime = primes[len(primes) - 1] + for i in range(last_prime, n): + if is_prime(i + 1): + primes.append(i + 1) + +def get_prime_factors(n : int) -> list: + factors = [] + init_primes_to(n) + for p in primes: + if p == n: + return [(p, 1)] + elif p > n: + break + exponent : int = 0 + while n % p == 0: + exponent += 1 + n /= p + if exponent > 0: + factors.append((p, exponent)) + return factors + +def get_fraction(f : float, mx_denom=100) -> tuple: + ''' +In order to get past floating point weirdness, we have to put the float +into a string first. + ''' + # assert(type(f) == Fraction) + # return f.as_integer_ratio() + return Fraction(f).limit_denominator(max_denominator=mx_denom).as_integer_ratio() + +def minimal_scaling_factor(vec : np.matrix) -> int: + ''' +This does not work if the data type of the matrix is `float` + ''' + primes = {} + for elem in vec: + n, d = get_fraction(elem[0, 0]) + factors = get_prime_factors(d) + for p, e in factors: + primes[p] = max(primes[p], e) if p in primes else e + scale_factor = 1 + for prime, exponent in primes.items(): + scale_factor *= prime ** exponent + return scale_factor + +def get_nullvectors(R : np.matrix) -> list: + ''' +Gets the nullvectors of matrix R (where R has all positive integers) +and ensures that all of the nullvectors are also of type int. + ''' + return null_space(R) # Todo: turn into list of columns + +def get_cycle_vectors(R : np.matrix, num=5): + ''' +Any positive integer linear combination of the nullvectors of R are the cycles +of the graph. This gives us a set of `num` *reasonably small* cycle vectors. + ''' + nv = get_nullvectors(R) + cycles = [v for v in nv] + n = len(nv) + # Add linear combinations + # TODO: add support for combinations beyond that + print("[WARNING] Wayfarer only supports \"first level\" cycle detection currently. This means only linear combinations of null vectors with coefficients equal to 1 or 0") + # m_fac = 1 + # We already have each individual null vector in `cycles` + for i in range(2, n): + for j in range(0, i): + # TODO: I think this works + if len(cycles) >= num: + return cycles + vsum = sum([v for v in nv[j::i]]) + cycles.append(vsum) + # while len(cycles) < num: + # pass + return cycles + +def get_ordered_reactions(crn : Crn) -> list: + ''' +Orders the reactions by rate constant, returns a list of their +indexes in the crn's `transitions` member list + ''' + # Yeah this is messy, but whatever + sortable_transitions = [SortableTransition(t, 0) for t in crn.transitions] + for i in range(len(sortable_transitions)): + sortable_transitions[i].index = i + sortable_transitions.sort(reverse=True) + return [st.index for st in sortable_transitions] + +def cycles_from_cycle_vectors(vecs : list, crn : Crn) -> list: + ''' +Creates cycles from cycle vectors + ''' + cycles = [] + sorted_transitions = get_ordered_reactions(crn) + for v in vecs: + # Rather than combinatorially expand to all possible + # just get those with the most probable transitions + # first, since they are most likely to be probable + v_counter = v.copy() + transitions = [] + while not np.all(v_counter == 0): + for idx in sorted_transitions: + transitions.append(crn.transitions[i]) + v_counter[idx] -= 1 + cycles.append(Cycle(transitions, crn)) + + return cycles diff --git a/src/dependency.py b/src/dependency.py index 072b57a..8b5f609 100644 --- a/src/dependency.py +++ b/src/dependency.py @@ -6,6 +6,13 @@ from subspace import * from util import * +# from stormpy import Rational + +def null(A : np.matrix, eps : float=1e-7): + u, s, vh = np.linalg.svd(A, full_matrices=1, compute_uv=1) + null_space = np.compress(s <= eps, vh, axis=0) + return null_space.H + class Node: ''' A Node in the dependency graph @@ -145,9 +152,9 @@ def __init__(self, ragtimer_lines, agnostic=False): self.graph_root = self.create_graph(change) self.create_reaction_levels() - def create_offset_vector(self, sn : Subspace, s0 : Subspace = None): + def create_offset_vector(self, sn : Subspace, s0 : Subspace | None = None): ''' - Creates an offset vector, f, s.t. f \in S0 and f minimizes the distance from + Creates an offset vector, f, s.t. $f \\in S0$ and f minimizes the distance from Sn to Ss (the solution space). The process is as follows: @@ -155,14 +162,17 @@ def create_offset_vector(self, sn : Subspace, s0 : Subspace = None): 2. Then, we must find the shortest distance from Sn to this intersection ''' sa = self.particular_solution - self.init_state + Mpsivecs = self.sat_basis.copy() + if len(Mpsivecs) == 0: + return sa + # print(Avecs) + Mpsi = np.column_stack(Mpsivecs) + # print(A) + if s0 is None: # If this is true we can short circuit knowing that the # shortest distance is contained in S0 already - Avecs = self.sat_basis.copy() - for t in sn.transitions: - Avecs.append(-np.matrix(t.vector).T) - A = np.column_stack(Avecs) - offset_nonprojected = sa - A * np.linalg.pinv(A.T * A, rcond=1e-3) * A.T * sa + offset_nonprojected = sa - Mpsi * np.linalg.pinv(Mpsi.T * Mpsi, rcond=1e-3) * Mpsi.T * sa # Because the pinv is calculated numerically, on some models where we should # get an offset of zero-vector, we get something like [1e-14, 1e-15, ...]. These # perturbations actually affect the search distance, so we will just zero it here. @@ -174,30 +184,22 @@ def create_offset_vector(self, sn : Subspace, s0 : Subspace = None): return zeros return offset_nonprojected else: - # In this case, we wish to solve the equation - # Snx + s0 = P0(Ssy + sp - s0) + s0 (project the solution space onto the OFFSET S0) - # Snx = P0(Ssy + sp - s0) - # [Sn -P0Ss]v = P0(sp - s0) - saP = s0.P * sa - Avecs = [] - for v in self.sat_basis: - Avecs.append(s0.P * v) - for t in sn.transitions: - Avecs.append(-np.matrix(t.vector).T) - A = np.column_stack(Avecs) - offset_nonprojected = saP - A * np.linalg.pinv(A.T * A, rcond=1e-3) * A.T * saP - # Because the pinv is calculated numerically, on some models where we should - # get an offset of zero-vector, we get something like [1e-14, 1e-15, ...]. These - # perturbations actually affect the search distance, so we will just zero it here. - # TODO: how to handle numerical innacuracies in the actual subspace construction? - # MAKE SURE TO NOTE THIS IN THE PAPER - zeros = np.zeros(offset_nonprojected.shape) - if np.isclose(offset_nonprojected, zeros).all(): - print("Info: offset is zero vector. This means that the smallest subspace already intersects the solution space with no offset.") - print("\tThis does NOT mean that the solution space is reachable ONLY FROM THOSE REACTIONS from the initial state.") + b = self.particular_solution - 2 * self.init_state + # print(s0.M, Mpsi) + A = np.block([s0.M, -Mpsi]) + # print(A) + xy, _, _, _ = np.linalg.lstsq(A, b, rcond=None) + x = xy[:s0.M.shape[1]] + # y = xy[s0.M.shape[1]:] + + f = s0.M * x + self.init_state + f = np.round(f) + print(f) + zeros = np.zeros(f.shape) + if np.isclose(f, zeros).all(): # This also short-circuits the next projection step return zeros - return offset_nonprojected + return f def create_graph(self, change, level = 0): ''' @@ -292,14 +294,18 @@ def create_reaction_levels(self): self.create_reaction_level(node) self.reaction_levels.reverse() + # Should do minimal root distance!!!! def create_reaction_level(self, node : Node): if node is None: return -1 if node.children is None or len(node.children) == 0: + # print(f"Node {node.reaction.name} has level 0") self.declare_reaction_at_level(node.reaction, 0) return 0 + # print(f"Node {node.reaction.name} has children {' '.join([c.reaction.name for c in node.children])}") successor_levels = [self.create_reaction_level(n) for n in node.children] level = 1 + min(successor_levels) + # print(f"Node {node.reaction.name} has level {level}") self.declare_reaction_at_level(node.reaction, level) return level @@ -310,7 +316,9 @@ def declare_reaction_at_level(self, reaction, level): ''' while len(self.reaction_levels) <= level: self.reaction_levels.append([]) - self.reaction_levels[level].append(reaction) + level_list = self.reaction_levels[level] + if reaction not in level_list: + level_list.append(reaction) def declare_producer(self, reaction, species): ''' @@ -353,6 +361,7 @@ def create_subspaces(self, crn): level = self.reaction_levels[lev_idx] for reaction in level: transition = crn.find_transition_by_name(reaction.name) + transition.in_s0 = True available_reactions.append(transition) indecies.append(len(available_reactions) - 1) # The last index is unnecessary as it is just the last available index diff --git a/src/distance.py b/src/distance.py index c9af66f..e2a222b 100644 --- a/src/distance.py +++ b/src/distance.py @@ -1,4 +1,5 @@ import numpy as np +# from stormpy import Rational from crn import * diff --git a/src/main.py b/src/main.py index 122753b..c67da09 100755 --- a/src/main.py +++ b/src/main.py @@ -11,8 +11,8 @@ store_traces = False -def basic_priority(filename, num): - crn = parse_ragtimer(filename) +def basic_priority(filename, bound_fname, num): + crn = parse_ragtimer(filename, bound_fname) print("========================================================") print("Targeted Exploration (just distance)") @@ -29,8 +29,8 @@ def basic_priority(filename, num): # end_time = time.time() # print(f"Total time {end_time - start_time} s") -def random(filename, num): - crn = parse_ragtimer(filename) +def random(filename, bound_fname, num): + crn = parse_ragtimer(filename, bound_fname) print("========================================================") print("Random Exploration") print("========================================================") @@ -39,7 +39,7 @@ def random(filename, num): end_time = time.time() print(f"Total time {end_time - start_time} s") -def subspace_priority(filename, num): +def subspace_priority(filename, bound_fname, num): dep, crn = parse_dependency_ragtimer(filename) print("========================================================") print("Targeted Exploration (Subspace)") @@ -49,8 +49,8 @@ def subspace_priority(filename, num): end_time = time.time() print(f"Total time {end_time - start_time} s") -def subspace_priority_solver(filename, num, time_bound, agnostic=False, piped=False, all_expand=False, single_order=False, use_rate_const=False): - dep, crn = parse_dependency_ragtimer(filename, agnostic=agnostic) +def subspace_priority_solver(filename, bound_fname, num, time_bound, agnostic=False, piped=False, all_expand=False, single_order=False, use_rate_const=False): + dep, crn = parse_dependency_ragtimer(filename, agnostic=agnostic, dimension_bounds_filename=bound_fname) print("========================================================") print("Targeted Exploration (Subspace - With Solver)") print("========================================================") @@ -66,10 +66,12 @@ def subspace_priority_solver(filename, num, time_bound, agnostic=False, piped=Fa if __name__=="__main__": parser = argparse.ArgumentParser( prog="wayfarer" - , description="wayfarer is a heuristic counterexample generator compatible with RAGTIMER files" + , description="wayfarer is a heuristic partial state space generator for lower bounds compatible with RAGTIMER files" , epilog="Developed at USU") parser.add_argument("-r", "--ragtimer", default=None, help="The name of the .ragtimer file to check") + parser.add_argument("-b", "--variable_bounds", default=None, + help="The name of the file that imposes upper bounds on the variables") parser.add_argument("-n", "--number", default=3, help="The number of either counterexamples or satisfying states to find, if without or with solver respectively") parser.add_argument("-p", "--primitive", action="store_true", @@ -97,23 +99,33 @@ def subspace_priority_solver(filename, num, time_bound, agnostic=False, piped=Fa help="When expanding states, expand ALL transitions rather than just those which the dependency graph specifies will get the shortest traces. This may help find higher probability traces.") parser.add_argument("-Q", "--rate_constant", action="store_true", help="Use rate constant in piped method.") + parser.add_argument("-F", "--rate_finder", default=None, + help="If reaction rate ought to be calculated by something other than the standard rate constant method, a Python file may be provided here in order to perform these custom calculations. Will look for a function called rate_finder(state : arraylike, rate_constant : float, reaction_name : str) -> float") + parser.add_argument("-U", "--upper", action="store_true", + help="Also compute upper bound") args = parser.parse_args() store_traces = args.traces if args.ragtimer is None: print("Missing args.") sys.exit(1) num = int(args.number) + + if args.rate_finder is not None: + parse_custom_rate_finder(args.rate_finder) + + SolverSettings.COMPUTE_UPPER_BOUND = args.upper if args.subspace: - subspace_priority(args.ragtimer, num) + subspace_priority(args.ragtimer, args.variable_bounds, num=num) if args.subspace_with_solver: t = None - if args.time is not None and not args.time.isnumeric(): + if args.time is not None and not args.time.replace(".", "").isnumeric(): print(f"Time bound {args.time} is invalid. Will ignore.") elif args.time is not None: - t = int(args.time) + t = float(args.time) subspace_priority_solver(args.ragtimer - , num + , args.variable_bounds + , num=num , time_bound=t , agnostic=args.agnostic , piped=args.piped @@ -122,12 +134,13 @@ def subspace_priority_solver(filename, num, time_bound, agnostic=False, piped=Fa if args.solver: t = None - if args.time is not None and not args.time.isnumeric(): + if args.time is not None and not args.time.replace(".", "").isnumeric(): print(f"Time bound {args.time} is invalid. Will ignore.") elif args.time is not None: - t = int(args.time) + t = float(args.time) subspace_priority_solver(args.ragtimer - , num + , args.variable_bounds + , num=num , time_bound=t , agnostic=args.agnostic , piped=args.piped @@ -136,8 +149,8 @@ def subspace_priority_solver(filename, num, time_bound, agnostic=False, piped=Fa , use_rate_const=args.rate_constant) if args.primitive: - basic_priority(args.ragtimer, num) + basic_priority(args.ragtimer, args.variable_bounds, num=num) if args.random: - random(args.ragtimer, num) + random(args.ragtimer, args.variable_bounds, num=num) diff --git a/src/parser.py b/src/parser.py index f79d782..678d24a 100644 --- a/src/parser.py +++ b/src/parser.py @@ -2,7 +2,32 @@ from distance import Bound, BoundTypes from dependency import * +# from stormpy import Rational + import sys +# Cursed +import importlib +import importlib.util + +# A dynamically loaded function pointer from a user file to customize how rates are found +custom_rate_finder = None + +def parse_custom_rate_finder(filename : str): + ''' +Loads `filename` (a python file) and looks for a function called rate_finder which is then +used to generate rates for a CRN + ''' + global custom_rate_finder # this should be the ONLY function which MODIFIES this + assert(filename.endswith(".py")) + module_name = filename.replace(".py", "").split("/").pop() + # mod = importlib.import_module(module_name) + module_spec = importlib.util.spec_from_file_location(module_name, filename) + mod = importlib.util.module_from_spec(module_spec) + module_spec.loader.exec_module(mod) + try: + custom_rate_finder = mod.rate_finder + except Exception as e: + print(f"[WARNING] Could not load custom rate finder. Reason: {e}") def check_valid_identifier(identifier): ''' @@ -28,7 +53,21 @@ def create_bound(bound_text): return Bound(0, BoundTypes.DONT_CARE) return Bound(bound_int, BoundTypes.EQUAL) -def create_transition(transition_line, species_idxes): +def get_val_quantity(s: str) -> tuple: + ''' +Returns a tuple (species_name, quantity) + ''' + split_line = s.split("*") + if len(split_line) == 1: + return split_line[0], 1 + elif len(split_line) == 2: + first, second = split_line + if first.isnumeric(): + return second, int(first) + else: + return first, int(second) + +def create_transition(transition_line, species_idxes, dimension_bounds=None): transition_info = transition_line.split("\t") sep_idx = transition_info.index(">") tname = transition_info[0] @@ -38,35 +77,47 @@ def create_transition(transition_line, species_idxes): transition_vector = np.array([0 for _ in range(len(species_idxes))]) always_enabled = False is_consumer = False - for reactant in reactants: + # Is 1.0 iff is reactant + rate_mul_vector = np.array([0.0 for _ in transition_vector]) + reactant_tuples = map(get_val_quantity, reactants) + for reactant, count in reactant_tuples: if reactant == "0": always_enabled = True break - transition_vector[species_idxes[reactant]] -= 1 - for product in products: + idx = species_idxes[reactant] + transition_vector[idx] -= count + rate_mul_vector[idx] = 1.0 + for product, count in map(get_val_quantity, products): if product == "0": is_consumer = True break transition_vector[species_idxes[product]] += 1 - # Is 1.0 iff is reactant - rate_mul_vector = np.array([float(elem < 0) for elem in transition_vector]) # The rate, from rate constant k and reactants A, B, is k * A^count(A) * B^count(B) - rate_finder = lambda state : rate_const * np.prod([state[i] ** rate_mul_vector[i] for i in range(len(rate_mul_vector))]) + rate_finder = None + if custom_rate_finder is None: + rate_finder = lambda state : rate_const * np.prod([state[i] ** rate_mul_vector[i] for i in range(len(rate_mul_vector))]) + else: + rate_finder = lambda state : custom_rate_finder(state, rate_const, tname) if always_enabled: return Transition(transition_vector , lambda state : True , rate_finder , tname - , rate_const) - reactant_idxes = [species_idxes[reactant] for reactant in reactants] + , rate_const + , dim_bounds=dimension_bounds) + reactant_idxes = [species_idxes[reactant] for reactant, _ in reactant_tuples] # Require all reactants to be strictly greater than zero return Transition(transition_vector , lambda state : np.all([state[i] > 0 for i in reactant_idxes]) , rate_finder , tname - , rate_const) + , rate_const + , dim_bounds=dimension_bounds) -def parse_ragtimer(filename): +def parse_ragtimer(filename, dimension_bounds_filename=None): + dimension_bounds = None if dimension_bounds_filename is None else create_dimension_bounds(dimension_bounds_filename) + # if dimension_bounds is not None: + # print(f"Got dimension bounds {dimension_bounds}") with open(filename, 'r') as rag: lines = rag.readlines() assert(len(lines) >= 4) @@ -77,14 +128,20 @@ def parse_ragtimer(filename): # Get the boundary condition (exact equal) bound_vals = lines[2].split("\t") boundary = [create_bound(bound_val) for bound_val in bound_vals] - transitions = [create_transition(line, species_idxes) for line in lines[3:]] + transitions = [create_transition(line, species_idxes, dimension_bounds) for line in lines[3:]] return Crn(transitions, boundary, init_state) -def parse_dependency_ragtimer(filename: str, agnostic : bool =False): +def create_dimension_bounds(dimension_bounds_filename : str) -> np.array: + with open(dimension_bounds_filename, 'r') as bfile: + lines = list(filter(lambda line: line != "", bfile.readlines())) + assert(len(lines) == 1) + return np.array([int(elem) for elem in lines[0].strip().split("\t")]) + +def parse_dependency_ragtimer(filename: str, agnostic : bool =False, dimension_bounds_filename : str | None = None): with open(filename, 'r') as rag: lines = rag.readlines() assert(len(lines) >= 4) - return DepGraph(lines, agnostic=agnostic), parse_ragtimer(filename) + return DepGraph(lines, agnostic=agnostic), parse_ragtimer(filename, dimension_bounds_filename=dimension_bounds_filename) def create_piped(crn : Crn, use_rc : bool = False): ''' diff --git a/src/rubound.py b/src/rubound.py index 8219752..d09c2f3 100644 --- a/src/rubound.py +++ b/src/rubound.py @@ -8,6 +8,7 @@ ''' import numpy as np +# from stormpy import Rational from crn import * diff --git a/src/solver.py b/src/solver.py index bbefe41..920f6ad 100644 --- a/src/solver.py +++ b/src/solver.py @@ -11,10 +11,11 @@ from stormpy import SparseMatrixBuilder, StateLabeling, SparseModelComponents import stormpy -DESIRED_NUMBER_STATES=2 -ABSORBING_INDEX=0 - -PRINT_FREQUENCY=100000 +class SolverSettings: + DESIRED_NUMBER_STATES=2 + ABSORBING_INDEX=0 + PRINT_FREQUENCY=100000 + COMPUTE_UPPER_BOUND=False all_states = [] state_ids = {} @@ -100,9 +101,9 @@ def assert_all_entries_correct(self): continue max_entry = max(self.from_list[i]) max_rate = max_entry.val - if not self.exit_rates[i] >= max_rate: + if self.exit_rates[i] is not None and not self.exit_rates[i] >= max_rate: print(f"Error: {self.exit_rates[i]} < {max_rate} (state index {i})") - assert(self.exit_rates[i] >= max_rate or math.isclose(max_rate, self.exit_rates[i])) + assert(self.exit_rates[i] is None or (self.exit_rates[i] >= max_rate or math.isclose(max_rate, self.exit_rates[i]))) def min_probability_subsp(crn, dep, number=1, print_when_done=False, write_when_done=False, time_bound=None, expand_all_states=False, single_order=False): global all_states @@ -112,7 +113,7 @@ def min_probability_subsp(crn, dep, number=1, print_when_done=False, write_when_ all_states = [] # Add the absorbing state matrixBuilder = RandomAccessSparseMatrixBuilder() - last_index = ABSORBING_INDEX + 1 + last_index = SolverSettings.ABSORBING_INDEX + 1 all_states.append(None) # Other stuff boundary = crn.boundary @@ -132,7 +133,7 @@ def min_probability_subsp(crn, dep, number=1, print_when_done=False, write_when_ num_explored = 0 while (not pq.empty()) and num_satstates < number: num_explored += 1 - if num_explored % PRINT_FREQUENCY == 0: + if num_explored % SolverSettings.PRINT_FREQUENCY == 0: print(f"Explored {num_explored} states. Have {num_satstates} satisfying") curr_state_data = pq.get() # print(f"Exploring state with index {curr_state_data.idx}") @@ -152,7 +153,8 @@ def min_probability_subsp(crn, dep, number=1, print_when_done=False, write_when_ # Total full rate: the total rate of all POSSIBLE enabled transitions from this state. successors, total_expanded_rate = curr_state_data.successors(all_successors=expand_all_states) total_full_rate = curr_state_data.get_total_outgoing_rate() - assert(total_full_rate >= total_expanded_rate) + # print(total_full_rate, total_expanded_rate) + assert(total_full_rate + 1e-5 >= total_expanded_rate) if len(successors) == 0: print("No successors") # Introduce a self-loop @@ -180,7 +182,10 @@ def min_probability_subsp(crn, dep, number=1, print_when_done=False, write_when_ pq.put(s) # If this state already exists, use the state data we already have else: + # update reachability for re-explored states + s_old = s s = all_states[state_ids[next_state_tuple]] + s.reach += s_old.reach assert(s.idx is not None) # Place the transition in the matrix matrixBuilder.add_next_value(curr_state_data.idx, s.idx, rate) @@ -216,12 +221,12 @@ def finalize_and_check(matrixBuilder : RandomAccessSparseMatrixBuilder, satisfyi if satisfies(state.vec, crn.boundary): # print(f"Found satisfying state {tuple(curr_state)}") num_perim_satstates += 1 - sat_states.append(curr_state_data.idx) + satisfying_state_idxs.append(state.idx) # We will create a self-loop later, so declare the total exit rate as 1.0 - matrixBuilder.add_exit_rate(curr_state_data.idx, 1.0) - matrixBuilder.add_next_value(curr_state_data.idx, curr_state_data.idx, 1.0) - deadlock_idxs.append(curr_state_data.idx) - curr_state_data.perimeter = False + matrixBuilder.add_exit_rate(state.idx, 1.0) + matrixBuilder.add_next_value(state.idx, state.idx, 1.0) + deadlock_idxs.append(state.idx) + state.perimeter = False continue # Expand the state and create transitions ONLY TO EXISTING STATES successors, total_exit_rate = state.successors(True) @@ -269,5 +274,12 @@ def finalize_and_check(matrixBuilder : RandomAccessSparseMatrixBuilder, satisfyi env = stormpy.Environment() env.solver_environment.native_solver_environment.precision = stormpy.Rational(1e-100) result = stormpy.check_model_sparse(model, prop, only_initial_states=True) - assert(result.min >= 0.0 and result.max <= 1.0) print(f"Pmin = {result.at(1)}") + assert(result.min + 1e-6 >= 0.0 and result.max <= 1.0 + 1e-6) + if SolverSettings.COMPUTE_UPPER_BOUND: + # Upper bound propert + chk_property_upper = f"P=? [ true U{prop_bound} \"satisfy\" | \"absorbing\" ]" + # print(f"Checking upper bound with property `{chk_property_upper}`") + prop_upper = stormpy.parse_properties(chk_property_upper)[0] + result_upper = stormpy.check_model_sparse(model, prop_upper, only_initial_states=True) + print(f"Pmax = {result_upper.at(1)}") diff --git a/src/subspace.py b/src/subspace.py index 776d55e..420d047 100644 --- a/src/subspace.py +++ b/src/subspace.py @@ -9,22 +9,27 @@ # if VERIFY: # from nagini_contracts.contracts import * # from typing import List +from fractions import Fraction from distance import vass_distance from crn import * +from util import to_frac_matrix + +# from stormpy import Rational DONT_CARE = -1 class Subspace: - mask = None # (target > DONT_CARE).astype(float) + mask = None # (target > DONT_CARE).astype(float) # This parillelipiped is the piped matrix containing the normalized # reaction vector scaled by the reaction rate, assuming a constant rate. # Then, piped_inv is the inverse of that Piped matrix, computed only once # for brevity and optimization. - piped = None # Stored for completeness + piped = None # Stored for completeness piped_inv = None # Assumes that piped ** -1 = piped_inv def initialize_piped(piped_matrix : np.matrix) -> None: Subspace.piped = piped_matrix + to_frac_matrix(Subspace.piped) Subspace.piped_inv = np.linalg.pinv(piped_matrix) # @Pure @@ -42,8 +47,8 @@ def norm(vec): # -> float: # Requires(len(vec) == len(Subspace.mask)) # Ensures(Result() >= 0.0) if Subspace.piped_inv is not None: - return float(np.linalg.norm(Subspace.piped_inv * vec)) - return float(np.linalg.norm(np.multiply(vec, Subspace.mask))) + return Fraction(np.linalg.norm(Subspace.piped_inv * vec, ord=1)).limit_denominator(300) + return Fraction(np.linalg.norm(np.multiply(vec, Subspace.mask), ord=1)).limit_denominator(300) # Type of elements in transitions: crn.Transition @@ -64,6 +69,7 @@ def __init__(self, transitions, excluded_transitions, last_layer = None): # print(basis_vectors[0]) # TODO: make sure we're appending to the right axis A = np.column_stack(basis_vectors) # np.matrix(basis_vectors[0].append(basis_vectors[1:], axis=1)) + to_frac_matrix(A) self.M = A # print(A) # Use the pseudoinverse since with rectangular matrices, @@ -89,7 +95,6 @@ def get_update_vectors(self, crn=None): # @Pure def contains(self, other, test_vec): # -> bool: # Check if contains - # Requires(type(State.init) == np.matrix) # Requires(len(test_vec) == len(Subspace.mask)) # Ensures(Implies(Result(), self.dist(test_vec) >= other.dist(test_vec))) if self.rank < other.rank: @@ -117,6 +122,11 @@ class State: init : np.matrix = None crn : Crn = None total_offset : np.matrix = None + # Cycle and commute stuff. There is a list of known useful cycles, orthocycles, + # and trivially commutable transitions + orthocycles : list = [] + non_orthocycles : list = [] + commutable_transitions : list = [] # @staticmethod def initialize_static_vars(crn, dep, single_order=False): if not single_order: @@ -127,6 +137,9 @@ def initialize_static_vars(crn, dep, single_order=False): State.init = np.matrix(crn.init_state).T State.target = np.matrix([b.to_num() for b in crn.boundary]).T Subspace.mask = np.matrix([b.to_mask() for b in crn.boundary]).T + to_frac_matrix(State.init) + to_frac_matrix(State.target) + to_frac_matrix(Subspace.mask) State.crn = crn if len(State.subspaces) == 0: State.total_offset = State.init @@ -136,9 +149,10 @@ def initialize_static_vars(crn, dep, single_order=False): # Result vector must be projected on s0.P else: State.total_offset = State.init + dep.create_offset_vector(State.subspaces[len(State.subspaces) - 1], State.subspaces[0]) + to_frac_matrix(State.total_offset) print(f"{dep}") - def __init__(self, vec, idx=None): + def __init__(self, vec, idx=None, reach=1.0): ''' Constructor for a new State element. Members within the State class: 1. vec (type: np.matrix) : the actual vector representing the state values @@ -150,43 +164,48 @@ def __init__(self, vec, idx=None): 4. epsilon (type list(int)) : The list of nonzero subspace distances, starting with the largest and working towards the smallest ''' - # Requires(type(State.init) == np.matrix) - # Requires(type(State.target) == np.matrix) - # Requires(type(State.subspaces) == List[Subspace]) # # Requires(Forall(int, lambda i : Implies(i > 0 and i < len(State.subspaces), State.subspaces[i].contains(State.subspaces[i - 1], State.init)))) # Ensures(self.order >= -1) # Ensures(len(self.epsilon) == len(State.subspaces) + 1) self.vec = vec self.vecm = np.matrix(vec).T + to_frac_matrix(self.vecm) self.adj = self.vecm - State.total_offset # State.init self.order : int = 0 self.__compute_order() self.perimeter = True self.idx = idx + self.sbsp = State.subspaces[0] if len(State.subspaces) > 0 else None + self.reach = reach def __compute_order(self): ''' Computes the order and epsilon vector of the state ''' - # Requires(type(State.init) == np.matrix) - # Requires(type(State.target) == np.matrix) - # Ensures(type(self.order) == int and self.order >= -1) - # Ensures(type(self.epsilon == list[float])) # Ensures(len(self.epsilon) >= 1) # Ensures(Forall(int, lambda i : (Implies(i > 0 and i < len(State.subspaces), self.epsilon[i] >= self.epsilon[i - 1])))) dist_to_target = Subspace.norm(self.vecm - State.target) - if dist_to_target == 0.0: - self.epsilon : list = [0.0] + if dist_to_target == 0: + self.epsilon : list = [0] self.order = -1 return self.epsilon = [dist_to_target] self.order = 0 - for s in State.subspaces: + # print(self.vec, end=": ") + for s in reversed(State.subspaces): + # print(s.rank, end=",") ep = s.dist(self.adj) - if ep == 0: + # For some reason the floating point thing has some issues + #if ep == 0: + #if ep > 1e-14 and ep < 1e-12: + # print(f"Gotcha! {ep} on state {self.vec} for {s}") + if ep < 1e-8: # To account for floating point error + self.sbsp = s + # print() return self.epsilon.insert(0, ep) self.order += 1 + # print() def get_total_outgoing_rate(self): ''' @@ -200,7 +219,11 @@ def get_total_outgoing_rate(self): total_rate = 0.0 for t in transitions: if t.enabled(vec): - total_rate += t.rate_finder(vec) + rate = t.rate_finder(vec) + assert(rate >= 0) + if rate == 0.0: + continue + total_rate += rate return total_rate def successors(self, only_tuples : bool = False, all_successors : bool = False): # -> tuple: @@ -208,12 +231,7 @@ def successors(self, only_tuples : bool = False, all_successors : bool = False): Only returns the successors using the vectors in the dependency graph that get us closer to the target. ''' - # Requires(type(State.init) == np.matrix) - # Requires(type(State.target) == np.matrix) - # Ensures(type(Result()) == tuple) # Ensures(len(Result()) == 2) - # Ensures(type(Result()[0]) == list) - # Ensures(type(Result()[1]) == float) # Ensures(Result()[1] >= 0.0) # If we get the successors, we are no longer a perimeter state @@ -226,17 +244,22 @@ def successors(self, only_tuples : bool = False, all_successors : bool = False): update_vectors = State.crn.transitions else: # TODO: why is this IndexError'ing on some models? - subspace = State.subspaces[max(0, len(State.subspaces) - (self.order + 2))] + subspace = self.sbsp # State.subspaces[max(0, len(State.subspaces) - (self.order + 2))] + # print(max(0, len(State.subspaces) - (self.order + 2))) if all_successors: # If the CRN variable is passed into get_update_vectors, all successors are returned - update_vectors = subspace.get_update_vectors(State.crn) + update_vectors = State.crn.transitions # subspace.get_update_vectors(State.crn) else: update_vectors = subspace.get_update_vectors() # State.crn) - # print(f"Update vectors {[str(vec) for vec in update_vectors]}") + # print(f"Update vectors {[str(vec.name) for vec in update_vectors]}") for t in update_vectors: # print(f"Update vector: {t.name} vec {t.vector}...", end="") if t.enabled(self.vec): rate = t.rate_finder(self.vec) + assert(rate >= 0.0) + # if we get a zero rate we can ignore things + if rate == 0.0: + continue total_outgoing_rate += rate if only_tuples: succ.append((tuple(self.vec + t.vector), rate)) @@ -263,24 +286,37 @@ def successors(self, only_tuples : bool = False, all_successors : bool = False): if t.enabled(self.vec): rate = t.rate_finder(self.vec) total_outgoing_rate += rate + # This assumes these states are new. State re-exploring is handled outside of + # this function (in the main exploration function) + if not only_tuples: + for next_state, rate in succ: + next_state.reach = self.reach * (rate / total_outgoing_rate) # print([s[0].vec for s in succ]) # print(f"For state {self.vec}, successors are {[state.vec for state, rate in succ]}") return succ, total_outgoing_rate + # TODO: I think there may be a bug here in what is being compared + # Comparators. ONLY COMPARES THE ORDER AND THE LOWEST VALUE FOR EPSILON # @Pure def __gt__(self, other): # Requires(len(self.epsilon) == len(other.epsilon)) # Requires(len(self.epsilon) > 0) # return self.epsilon[len(self.epsilon) - 1] > other.epsilon[len(other.epsilon) - 1] - return self.order > other.order or (self.order == other.order and self.epsilon[0] > other.epsilon[0]) + return self.order > other.order or \ + (self.order == other.order and self.epsilon[0] > other.epsilon[0]) or \ + (self.order == other.order and self.epsilon[0] == other.epsilon[0] and \ + self.reach < other.reach) # @Pure def __lt__(self, other): # Requires(len(self.epsilon) == len(other.epsilon)) # Requires(len(self.epsilon) > 0) # return self.epsilon[len(self.epsilon) - 1] < other.epsilon[len(other.epsilon) - 1] - return self.order < other.order or (self.order == other.order and self.epsilon[0] < other.epsilon[0]) + return self.order < other.order or \ + (self.order == other.order and self.epsilon[0] < other.epsilon[0]) or \ + (self.order == other.order and self.epsilon[0] == other.epsilon[0] and \ + self.reach > other.reach) # @Pure def __le__(self, other): @@ -304,7 +340,7 @@ def __eq__(self, other): # Requires(len(self.epsilon) == len(other.epsilon)) # Requires(len(self.epsilon) > 0) # return self.epsilon[len(self.epsilon) - 1] == other.epsilon[len(other.epsilon) - 1] - return self.order == other.order and self.epsilon[0] == other.epsilon[0] + return self.order == other.order and self.epsilon[0] == other.epsilon[0] and self.reach == other.reach # Strong equality means we are the same state # @Pure diff --git a/src/util.py b/src/util.py index 649ad76..e591c97 100644 --- a/src/util.py +++ b/src/util.py @@ -1,44 +1,55 @@ import numpy as np from numpy.linalg import svd +from fractions import Fraction + # Taken from https://scipy-cookbook.readthedocs.io/items/RankNullspace.html def nullspace(A, atol=1e-13, rtol=0): - """Compute an approximate basis for the nullspace of A. - - The algorithm used by this function is based on the singular value - decomposition of `A`. - - Parameters - ---------- - A : ndarray - A should be at most 2-D. A 1-D array with length k will be treated - as a 2-D with shape (1, k) - atol : float - The absolute tolerance for a zero singular value. Singular values - smaller than `atol` are considered to be zero. - rtol : float - The relative tolerance. Singular values less than rtol*smax are - considered to be zero, where smax is the largest singular value. - - If both `atol` and `rtol` are positive, the combined tolerance is the - maximum of the two; that is:: - tol = max(atol, rtol * smax) - Singular values smaller than `tol` are considered to be zero. - - Return value - ------------ - ns : ndarray - If `A` is an array with shape (m, k), then `ns` will be an array - with shape (k, n), where n is the estimated dimension of the - nullspace of `A`. The columns of `ns` are a basis for the - nullspace; each element in numpy.dot(A, ns) will be approximately - zero. - """ - - A = np.atleast_2d(A) - u, s, vh = svd(A) - tol = max(atol, rtol * s[0]) - nnz = (s >= tol).sum() - ns = vh[nnz:].conj().T - return ns + """Compute an approximate basis for the nullspace of A. + + The algorithm used by this function is based on the singular value + decomposition of `A`. + + Parameters + ---------- + A : ndarray + A should be at most 2-D. A 1-D array with length k will be treated + as a 2-D with shape (1, k) + atol : float + The absolute tolerance for a zero singular value. Singular values + smaller than `atol` are considered to be zero. + rtol : float + The relative tolerance. Singular values less than rtol*smax are + considered to be zero, where smax is the largest singular value. + + If both `atol` and `rtol` are positive, the combined tolerance is the + maximum of the two; that is:: + tol = max(atol, rtol * smax) + Singular values smaller than `tol` are considered to be zero. + + Return value + ------------ + ns : ndarray + If `A` is an array with shape (m, k), then `ns` will be an array + with shape (k, n), where n is the estimated dimension of the + nullspace of `A`. The columns of `ns` are a basis for the + nullspace; each element in numpy.dot(A, ns) will be approximately + zero. + """ + + A = np.atleast_2d(A) + _, s, vh = svd(A) + tol = max(atol, rtol * s[0]) + nnz = (s >= tol).sum() + ns = vh[nnz:].conj().T + return ns + +def to_frac_matrix(m : np.matrix): + ''' +Creates a matrix where all elements are Fraction types + ''' + I, J = m.shape + for i in range(I): + for j in range(J): + m[i, j] = Fraction(m[i, j]).limit_denominator(100) diff --git a/test.py b/test.py new file mode 100755 index 0000000..69ace71 --- /dev/null +++ b/test.py @@ -0,0 +1,77 @@ +#!/usr/bin/env python3 +import time +import subprocess +import re +import sys +import os + +models = [ + ("models/2react/2react", 100) + , ("models/6react/6react", 100) + , ("models/8react/8react", 20) + , ("models/12react/12react", 10) +] + +def getResultFromOutput(output : str) -> float: + pMin : float = 0.0 + for line in output.split("\n"): + if line.startswith("Pmin = "): + pMin = float(line.replace("Pmin = ", "")) + break + return pMin + +def getFoundStatesFromOutput(output : str) -> int: + foundStates : int = 0 + matches = re.findall(r'.+Found ([0-9]+) satisfying states.', output) + print(len(matches)) + print(output) + assert(len(matches) == 1) + return matches[0] + +def runTest(fileBaseName : str, timeBound : int, num : int = 2, agnostic : bool = False, piped : bool = False, sop : bool = False) -> None: + print(f"Testing Wayfarer with file {fileBaseName} searching for {num} sat states") + modelFile = f"{fileBaseName}.sm" + ragtimerFile = f"{fileBaseName}.ragtimer" + propFile = f"{fileBaseName}.csl" + # Test wayfarer + startTime = time.time() + cmdArray = None + if sop: + print("Testing single order priority") + cmdArray = ["./wayfarer", "-r", ragtimerFile, "-V", "-T", f"{timeBound}", "-n", f"{num}"] + else: + cmdArray = ["./wayfarer", "-r", ragtimerFile, "-S", "-T", f"{timeBound}", "-n", f"{num}"] + if agnostic: + cmdArray.append("--agnostic") + if piped: + cmdArray.append("--piped") + + p = subprocess.run(cmdArray, capture_output=True) + endTime = time.time() + results = p.stdout.decode("utf-8") + return getFoundStatesFromOutput(results), getResultFromOutput(results), endTime - startTime + +if __name__=="__main__": + postfix="/" + if "WAYFARER_TEST_POSTFIX" in os.environ: + postfix = f"{os.environ['WAYFARER_TEST_POSTFIX']}/" + print(f"Putting results in results/{postfix}") + agnostic = "-a" in sys.argv or "--agnostic" in sys.argv + sop = "-s" in sys.argv or "--sop" in sys.argv + if agnostic: + print("Initial state agnostic") + if sop: + print("Single order priority") + piped = "--piped" in sys.argv + if piped: + print("Piped method") + for baseName, tBound in models: + with open(f"results/{postfix}{baseName.split('/').pop()}.csv", "w") as f: + f.write("Number Sat Requested,Number Sat Found,Probability Min,Duration\n") + for p in range(14): + num = 2 ** p + #for p in range(4): + # num = 10 ** p + numFound, pMin, duration = runTest(baseName, tBound, num=num, agnostic=agnostic, piped=piped, sop=sop) + print(f"Took {duration} s and found pMin={pMin} ({numFound} actual found satisfying states)") + f.write(f"{num},{numFound},{pMin},{duration}\n") diff --git a/wayfarer b/wayfarer new file mode 120000 index 0000000..a8c92b0 --- /dev/null +++ b/wayfarer @@ -0,0 +1 @@ +src/main.py \ No newline at end of file