diff --git a/main_package/classes/amalgamated_cims.py b/main_package/classes/amalgamated_cims.py index b02b308..470cd82 100644 --- a/main_package/classes/amalgamated_cims.py +++ b/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) diff --git a/main_package/classes/conditional_intensity_matrix.py b/main_package/classes/conditional_intensity_matrix.py index 49ad3f0..82445dd 100644 --- a/main_package/classes/conditional_intensity_matrix.py +++ b/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 diff --git a/main_package/classes/network_graph.py b/main_package/classes/network_graph.py index 0b17fff..aa5aa0a 100644 --- a/main_package/classes/network_graph.py +++ b/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] diff --git a/main_package/classes/parameters_estimator.py b/main_package/classes/parameters_estimator.py index 5ba476c..4fbf03c 100644 --- a/main_package/classes/parameters_estimator.py +++ b/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) + diff --git a/main_package/classes/set_of_cims.py b/main_package/classes/set_of_cims.py index 547b9e1..e14d340 100644 --- a/main_package/classes/set_of_cims.py +++ b/main_package/classes/set_of_cims.py @@ -17,14 +17,33 @@ 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} - literal_index = "" - for node in self.ordered_parent_set: - literal_index = literal_index + str(dict_of_indexes[node]) - return int(literal_index, self.value) + #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)