Add Parameters EstimationAlgorithm

main_package/classes/amalgamated_cims.py

@@ -23,5 +23,8 @@ class AmalgamatedCims:
     def get_vars_order(self, node):
         return self.actual_cims[node][1]
 
-    def update_state_transition_for_matrix(self, node, dict_of_indxs, element_indx):
-        self.actual_cims[node]
+    def update_state_transition_for_matrix(self, node, dict_of_nodes_values, element_indx):
+        self.sets_of_cims[node].update_state_transition(dict_of_nodes_values, element_indx)
+
+    def update_state_residence_time_for_matrix(self, which_node, which_matrix, which_element, time):
+        self.sets_of_cims[which_node].update_state_residence_time(which_matrix, which_element, time)

main_package/classes/conditional_intensity_matrix.py

@@ -4,10 +4,25 @@ import numpy as np
 class ConditionalIntensityMatrix:
 
     def __init__(self, dimension):
-        self.state_residence_times = np.zeros(shape=(1, dimension))
+        self.state_residence_times = np.zeros(shape=dimension)
         self.state_transition_matrix = np.zeros(shape=(dimension, dimension), dtype=int)
         self.cim = np.zeros(shape=(dimension, dimension), dtype=float)
 
-    def update_state_transition_count(self, positions_list):
-        self.state_transition_matrix[positions_list[0]][positions_list[1]] = self.state_transition_matrix[positions_list[0]][positions_list[1]] + 1
+    def update_state_transition_count(self, element_indx):
+        #print(element_indx)
+        self.state_transition_matrix[element_indx[0]][element_indx[1]] = \
+            self.state_transition_matrix[element_indx[0]][element_indx[1]] + 1
+
+    def update_state_residence_time_for_state(self, state, time):
+        #print("Time updating In state", state, time)
+        self.state_residence_times[state] = self.state_residence_times[state] + time
+
+    def compute_cim_coefficients(self):
+        for i, row in enumerate(self.state_transition_matrix):
+            row_sum = 0.0
+            for j, elem in enumerate(row):
+                rate_coefficient = elem / self.state_residence_times[i]
+                self.cim[i][j] = rate_coefficient
+                row_sum = row_sum + rate_coefficient
+            self.cim[i][i] = -1 * row_sum
 

main_package/classes/network_graph.py

@@ -19,18 +19,13 @@ class NetworkGraph():
 
 
     def init_graph(self):
-        #self.sample_path.build_trajectories()
-        #self.sample_path.build_structure()
         self.add_nodes(self.graph_struct.list_of_nodes())
         self.add_edges(self.graph_struct.list_of_edges())
 
     def add_nodes(self, list_of_nodes):
         for indx, id in enumerate(list_of_nodes):
-            #print(indx, id)
             self.graph.add_node(id)
             nx.set_node_attributes(self.graph, {id:indx}, 'indx')
-        #for node in list(self.graph.nodes):
-            #print(node)
 
     def add_edges(self, list_of_edges):
         self.graph.add_edges_from(list_of_edges)
@@ -59,8 +54,11 @@ class NetworkGraph():
     def get_states_number(self):
         return self.graph_struct.get_states_number()
 
-    def get_node_by_index(self, node_id):
-        return self.graph_struct.get_node_indx(node_id)
+    def get_node_by_index(self, node_indx):
+        return self.graph_struct.get_node_id(node_indx)
+
+    def get_node_indx(self, node_id):
+        return nx.get_node_attributes(self.graph, 'indx')[node_id]
 
 
     

main_package/classes/parameters_estimator.py

@@ -1,6 +1,8 @@
+import os
+import time as tm
+
 import network_graph as ng
 import sample_path as sp
-import os
 import amalgamated_cims as acims
 
 
@@ -16,6 +18,62 @@ class ParametersEstimator:
                                                      self.net_graph.get_nodes(),
                                                              self.net_graph.get_ord_set_of_par_of_all_nodes())
 
+    def parameters_estimation(self):
+        print("Starting computing")
+        t0 = tm.time()
+        for indx, trajectory in enumerate(self.sample_path.trajectories):
+            self.parameters_estimation_single_trajectory(trajectory.get_trajectory())
+            #print("Finished Trajectory number", indx)
+        t1 = tm.time() - t0
+        print("Elapsed Time ", t1)
+
+    def parameters_estimation_single_trajectory(self, trajectory):
+        #print(type(trajectory[0][0]))
+        for indx, row in enumerate(trajectory):
+            if trajectory[indx][1] == -1:
+                break
+            if trajectory[indx + 1][1] != -1:
+                transition = self.find_transition(trajectory[indx], trajectory[indx + 1])
+                which_node = self.net_graph.get_node_by_index(transition[0])
+            # print(which_node)
+                which_matrix = self.which_matrix_to_update(row, which_node)
+                which_element = transition[1]
+                self.amalgamated_cims_struct.update_state_transition_for_matrix(which_node, which_matrix, which_element)
+
+            #changed_node = which_node
+            time = self.compute_time_delta(trajectory[indx], trajectory[indx + 1])
+
+            for node in self.net_graph.get_nodes():
+                #if node != changed_node:
+                    # print(node)
+                    which_node = node
+                    which_matrix = self.which_matrix_to_update(row, which_node)
+                    which_element = row[self.net_graph.get_node_indx(node) + 1]
+                    # print("State res time element " + str(which_element) + node)
+                    # print("State res time matrix indx" + str(which_matrix))
+                    self.amalgamated_cims_struct.update_state_residence_time_for_matrix(which_node, which_matrix, which_element, time)
+
+
+    def find_transition(self, current_row, next_row):
+        for indx in range(1, len(current_row)):
+            if current_row[indx] != next_row[indx]:
+                return [indx - 1, (current_row[indx], next_row[indx])]
+
+
+    def compute_time_delta(self, current_row, next_row):
+        return next_row[0] - current_row[0]
+
+    def which_matrix_to_update(self, current_row, node_id): # produce strutture {'X':1, 'Y':2} dove X e Y sono i parent di node_id
+        result = {}
+        parent_list = self.net_graph.get_parents_by_id(node_id)
+        for node in parent_list:
+            result[node] = current_row[self.net_graph.get_node_indx(node) + 1]
+        # print(result)
+        return result
+
+
+
+
 # Simple Test #
 os.getcwd()
 os.chdir('..')
@@ -32,3 +90,11 @@ pe = ParametersEstimator(s1, g1)
 pe.init_amalgamated_cims_struct()
 print(pe.amalgamated_cims_struct.get_set_of_cims('X').get_cims_number())
 print(pe.amalgamated_cims_struct.get_set_of_cims('Y').get_cims_number())
+#pe.parameters_estimation_single_trajectory(pe.sample_path.trajectories[0].get_trajectory())
+pe.parameters_estimation()
+for matrix in pe.amalgamated_cims_struct.get_set_of_cims('Y').actual_cims:
+    print(matrix.state_residence_times)
+    print(matrix.state_transition_matrix)
+    matrix.compute_cim_coefficients()
+    print(matrix.cim)
+

main_package/classes/set_of_cims.py

@@ -17,13 +17,32 @@ class SetOfCims:
         for indx, matrix in enumerate(self.actual_cims):
             self.actual_cims[indx] = cim.ConditionalIntensityMatrix(self.value)
 
+    def update_state_transition(self, dict_of_indexes, element_indx_tuple):
+        matrix_indx = self.indexes_converter(dict_of_indexes)
+        #print("Converted Indx SRT")
+        #print(matrix_indx)
+        self.actual_cims[matrix_indx].update_state_transition_count(element_indx_tuple)
+
+    def update_state_residence_time(self, which_matrix, which_element, time):
+        matrix_indx = self.indexes_converter(which_matrix)
+        #print("Converted Indx")
+        #print(matrix_indx)
+        #print("Updating Time for variable ",self.node_id)
+        self.actual_cims[matrix_indx].update_state_residence_time_for_state(which_element, time)
+
+
     def get_cims_number(self):
         return len(self.actual_cims)
 
     def indexes_converter(self, dict_of_indexes): # Si aspetta oggetti del tipo {X:1, Y:1, Z:0}
+        #print(dict_of_indexes)
+        if not dict_of_indexes:
+            return 0
+        else:
             literal_index = ""
             for node in self.ordered_parent_set:
                 literal_index = literal_index + str(dict_of_indexes[node])
+                #print(literal_index)
             return int(literal_index, self.value)