Redo how offset vector $\vec{f}$ is calculated

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

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

dependencies.txt

@@ -0,0 +1,2 @@
+numpy
+stormpy

src/crn.py

@@ -23,21 +23,50 @@ 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
 		self.rate_finder = rate_finder
 		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

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