Extended time sensitivity

neet/boolean/network.py

@@ -9,6 +9,7 @@
 from neet.python import long
 from .sensitivity import SensitivityMixin
 import copy
+import itertools as itt
 
 
 class BooleanNetwork(SensitivityMixin, UniformNetwork):
@@ -159,7 +160,7 @@ def subspace(self, indices, state=None):
                 else:
                     i += 1
 
-    def hamming_neighbors(self, state):
+    def hamming_neighbors(self, state, c=1):
         """
         Get all states that one unit of Hamming distance from a given state.
 
@@ -176,13 +177,38 @@ def hamming_neighbors(self, state):
         :type state: list, numpy.ndarray
         :return: a list of neighbors of the given state
         :raises ValueError: if the state is not in the network's state space
+        """
+        state_copy = copy.copy(state)
+
+        if state not in self:
+            raise ValueError('state is not in state space')
+
+        I_comb_iter = itt.combinations(range(self.size), c)
+
+        def c_hamming_neighbors(self, state, c):
+            try:
+                nxt = next(I_comb_iter)
+                XORed = copy.copy(state_copy)
+                for i in nxt:
+                    XORed[i] ^= 1
+                return XORed
+            except StopIteration:
+                return None
+
+        neighbors = list()
+        neighbor = c_hamming_neighbors(self, state, c)
+        while neighbor is not None:
+            neighbors.append(copy.copy(neighbor))
+            neighbor = c_hamming_neighbors(self, state, c)
+
         """
         if state not in self:
             raise ValueError('state is not in state space')
         neighbors = [None] * self.size
         for i in range(self.size):
             neighbors[i] = copy.copy(state)
             neighbors[i][i] ^= 1
+        """
         return neighbors
 
     def distance(self, a, b):

neet/boolean/sensitivity.py

@@ -8,6 +8,13 @@
 import copy
 import numpy as np
 import numpy.linalg as linalg
+import math
+import itertools as itt
+import matplotlib.pyplot as plt
+
+"""
+.. import matplotlib.pyplot as plt
+"""
 
 
 class SensitivityMixin(object):
@@ -35,7 +42,7 @@ class SensitivityMixin(object):
     network models.
     """
 
-    def sensitivity(self, state, transitions=None):
+    def sensitivity(self, state, transitions=None, timesteps=1):
         """
         Compute the Boolean sensitivity at a given network state.
 
@@ -83,17 +90,21 @@ def sensitivity(self, state, transitions=None):
         """
         encoder = self._unsafe_encode
         distance = self.distance
-        neighbors = self.hamming_neighbors(state)
-
-        nextState = self.update(state)
+        neighbors = self.hamming_neighbors(state, c=1)
+        #neighbors_copy = [neighbor.copy() for neighbor in neighbors]
+        for t in range(timesteps):
+            nextState = self.update(state)
 
         # count sum of differences found in neighbors of the original
         s = 0.
+
         for neighbor in neighbors:
-            if transitions is not None:
-                newState = transitions[encoder(neighbor)]
-            else:
-                newState = self._unsafe_update(neighbor)
+            for t in range(timesteps):
+                if transitions is not None:
+                    newState = transitions[encoder(neighbor)]
+                    neighbor = newState
+                else:
+                    newState = self._unsafe_update(neighbor)
             s += distance(newState, nextState)
 
         return s / self.size
@@ -438,7 +449,7 @@ def lambdaQ(self, **kwargs):
         Q = self.average_difference_matrix(**kwargs)
         return max(abs(linalg.eigvals(Q)))
 
-    def average_sensitivity(self, states=None, weights=None, calc_trans=True):
+    def average_sensitivity(self, states=None, weights=None, calc_trans=True, timesteps=1):
         """
         Calculate average Boolean network sensitivity, as defined in
         [Shmulevich2004]_.
@@ -479,8 +490,243 @@ def average_sensitivity(self, states=None, weights=None, calc_trans=True):
 
         .. seealso:: :func:`sensitivity`
         """
+        if(timesteps==1):
+            Q = self.average_difference_matrix(states=states, weights=weights, calc_trans=calc_trans)
+            return np.sum(Q) / self.size
+        else:
+            total = 0
+            if calc_trans:
+                decoder = self.decode
+                trans = list(map(decoder, self.transitions))
+            else:
+                trans = None
+            if states is not None:
+                # print("states: ", states)
+                for state in states:
+                    # print("state: ", state)
+                    val = self.sensitivity(state, trans, timesteps)
+                    # print(" val= ", val)
+                    total += val
+                avg = total/ (len(states))
+                return avg
+            else:
+                for state in self:
+                    # print("state: ", state)
+                    # for t in range(1, timesteps) :
+                    val = self.sensitivity(state, trans, timesteps)
+                    # print(" val= ", val)
+                    total += val
+                avg = total/ self.volume
+                return avg
+
+    def c_sensitivity(self, state, transitions=None, c=1):
+
+        assert (isinstance(c, int)),"c needs to be an integer"
+        assert(c >= 0), "the value of c needs to be greater than or equal to zero"
+        assert (c <= self.size),"the value of c needs to be between 0 and the size of the network"
+
+
+
+        """
+        C-Sensitivity modification of the regular sensitivity function.
+
+        The c-sensitivity of :math:`f(x_1, //ldots, x_n)` at :math:`x` is defined as the number of 
+        c-Hamming neighbors of :math:`x` on which the function value is different from its value on :math:`x`. That is,
+
+        :param state: a single network state
+        :type state: list, numpy.ndarray
+        :param transitions: precomputed state transitions (*optional*)
+        :type transitions: list, numpy.ndarray, None
+        :return: the C-sensitivity at the provided state
+
+        """
+
+        encoder = self._unsafe_encode
+        distance = self.distance
+        state_copy = copy.copy(state)
+        nextState = self.update(state)
+
+        """
+        Returns an iterator for each vector I which is a strict subset of {1,...,n} and where |I| = c
+        note: if c = 0, this will return an empty tuple
+        """
+
+        I_comb_iter = itt.combinations(range(self.size), c)
+
+        """ 
+        Generator function which returns a new hamming neighbor
+        Each hamming neighbor is simply the product of self.state XOR I
+        """
+        def c_hamming_neighbors(self, state, c):
+            try:
+                nxt = next(I_comb_iter)
+                XORed = copy.copy(state_copy)
+                for i in nxt:
+                    XORed[i] ^= 1
+                return XORed
+            except StopIteration:
+                return None
+
+
+        """ 
+        #OK, so I messed with the function and it's only ~kinda~ a generator function now...
+        It's automatically advanced in the for loop, which
+        acts as a "try: next(neighbors); catch StopIteration:". This behavior is built 
+        into Python and is idiomatic.
+        """ 
+
+        s = 0.
+        neighbors_copy = []#uses quite a bit of memory and should be removed from code once unit testing is completed
+        copy_counter = 0
+
+        neighbor = c_hamming_neighbors(self,state,c)
+        while neighbor is not None:
+            neighbors_copy.append(copy.copy(neighbor))
+            if transitions is not None:
+                newState = transitions[encoder(neighbor)]
+            else:
+                newState = self._unsafe_update(neighbor)
+
+            # the paper which describes c-sensitivity uses an indicator function
+            # instead of a distance function. Distance will be used here instead of the indicator.
+
+            if distance(newState, nextState) > 0:
+                if c == 0:
+                    # The distance between 0-hamming neighbors should
+                    # be 0, since a 0-hamming neighbor is just the
+                    # original state/vector
+                    print("This shouldn't print. 0-hamming neighbors should not diverge.")
+                    print("c: ", c)
+                    print("1. neighbor:  ", neighbors_copy[copy_counter])
+                    print("2. state:     ", state_copy,"\n")
+                    print("1. newState:  ", newState)
+                    print("2. nextState: ", nextState,"\n\n")
+                s += distance(newState, nextState)
+            copy_counter += 1
+            neighbor = c_hamming_neighbors(self,state,c)
+        return s / copy_counter
+
+
+    def average_c_sensitivity(self, states=None, calc_trans=True, c=1):
+
+        """
+        Simple acts as a for-loop which does some precomputation before generating 
+        all possible states of the network (maintaining topology and connections, 
+        just changing the initial node values to all possible combinations of active
+        nodes)
+        Each generated state's c-sensitivity is summed and then divided by the total
+        number of generated states.
+        """
+
+
+        """
+
+        :param states: a set of network states
+        :type states: list, numpy.ndarray
+        :param transitions: precomputed state transitions (*optional*)
+        :type transitions: list, numpy.ndarray, None
+        :param calc_trans: pre-compute all state transitions; ignored if
+                           ``states`` or ``weights`` is ``None``
+        :type calc_trans: bool
+        :return: the sensitivity averaged over all possible 
+                            states of the network
+
+        if states is not None:
+            num_states = 0
+            states = list(states)
+            #print("type of states: ", type(states))
+            #print("len of states: ", type(states))
+            #print("len of states: ", len(states))
+            #print("len of states: ", len(states))
+            if calc_trans:
+                decoder = self.decode
+                trans = list(map(decoder, self.transitions))
+            else:
+                trans = None
+
+        """
+
+
+        s = 0
+
+        # optionally pre-calculate transitions
+        if calc_trans:
+            decoder = self.decode
+            trans = list(map(decoder, self.transitions))
+        else:
+            trans = None
+
+
+        if states is not None:
+            # Iterate through all provided states and sum 
+            # the distances between them and their c-hamming-neighbors
+            # after a single time-step (one synchronous transition)
+            for state in states:
+                s += self.c_sensitivity(state, trans, c)
+
+            # get the average of s by diving the sum by the
+            # number of states considered
+            s = s / len(states)
+
+
+        else:
+            # Generate all possible states 
+            # and sum the distances between them and each of their
+            # c-hamming-neighbors in the next time step (after one 
+            # synchronous transition)
+            for state in self:
+                s += self.c_sensitivity(state, trans, c)
+
+
+            # Deprecated implementation I'm keeping around while we implement unit testing. Will be removed after
+            # equivalence/correctness of the above (2 line) implementation is established
+            """for n in range(self.size):
+                state_gen = itt.combinations(range(self.size),n)
+                for state in state_gen:
+                    state_array = [0 for x in range(self.size)]
+                    for index in state:
+                        state_array[index] = 1
+                    s += self.c_sensitivity(state_array, trans, c)"""
+
+            # Average s by diving the sum by the 
+            # total number of possible states (2^n)
+            # s = s / np.power(2, self.size)
+            s = s / self.volume
+            """ 
+            s is now the average C-Sensitivity of f and must lie in the interval [0, (n choose c)] 
+            where n is the size of the network.
+            """
+
+            upper_bound = math.factorial(self.size) / (math.factorial(c) * math.factorial(self.size - c))
+            if s > upper_bound or s < 0:
+                raise RuntimeError('This value of S should not be possible and the code is therefore wrong')
+
+        return s
+        #yield s / upper_bound # yields the normalized average c-sensitivity
+
+    def derrida_plot(self, max_c=None):#, transitions=None): #X-Axis = c value, Y-Axis = output of Average_c_sensitivity
+        if max_c is None:
+            max_c = self.size
+
+        plt.title('Derrida Plot')
+        y_vals = []
+
+        """if transitions is not None:
+            for x in range(max_c):
+                y_vals.append(self.Average_c_sensitivity(self, transitions))
+        else:
+            for x in range(max_c):
+                y_vals.append(self.Average_c_sensitivity(self))"""
+        for x in range(max_c):
+                y_vals.append(self.Average_c_sensitivity(states=None, calc_trans=True, c=x))            
+
+        print(y_vals)
+    def Extended_Time_Plot(self, max_timesteps=4, transitions=None): #X-Axis = c value, Y-Axis = output of Extended_Time
+        if max_c is None:
+            max_c = self.size
 
-        Q = self.average_difference_matrix(states=states, weights=weights,
-                                           calc_trans=calc_trans)
+        plt.title('Extended_Time_Plot')
+        y_vals = []
 
-        return np.sum(Q) / self.size
+        for x in range(max_timesteps):
+            y_vals.append(self.average_sensitivity(states=None, weights=None, calc_trans=True, timesteps=x))

neet/landscape.py

@@ -142,7 +142,7 @@ def landscape(self, index=None, pin=None, values=None):
 
         This function implicitly calls :attr:`clear_landscape`, so make sure to
         create a reference to :attr:`landscape_data` if landscape information
-        has previously been compute and you wish to keep it around.
+        has previously been computed and you wish to keep it around.
 
         .. rubric:: Basic Usage