diff --git a/cdmir/discovery/PBSCM/PB_SCM.py b/cdmir/discovery/constraint/PBSCM/PB_SCM.py similarity index 100% rename from cdmir/discovery/PBSCM/PB_SCM.py rename to cdmir/discovery/constraint/PBSCM/PB_SCM.py diff --git a/cdmir/discovery/PBSCM/__init__.py b/cdmir/discovery/constraint/PBSCM/__init__.py similarity index 100% rename from cdmir/discovery/PBSCM/__init__.py rename to cdmir/discovery/constraint/PBSCM/__init__.py diff --git a/cdmir/discovery/PBSCM/util.py b/cdmir/discovery/constraint/PBSCM/util.py similarity index 100% rename from cdmir/discovery/PBSCM/util.py rename to cdmir/discovery/constraint/PBSCM/util.py diff --git a/cdmir/effect/PBSCM_PGF/CCARankTest.py b/cdmir/discovery/constraint/PBSCM_PGF/CCARankTest.py similarity index 100% rename from cdmir/effect/PBSCM_PGF/CCARankTest.py rename to cdmir/discovery/constraint/PBSCM_PGF/CCARankTest.py diff --git a/cdmir/effect/PBSCM_PGF/PB_SCM_PGF.py b/cdmir/discovery/constraint/PBSCM_PGF/PB_SCM_PGF.py similarity index 98% rename from cdmir/effect/PBSCM_PGF/PB_SCM_PGF.py rename to cdmir/discovery/constraint/PBSCM_PGF/PB_SCM_PGF.py index 3ed7bdc..0fd431f 100644 --- a/cdmir/effect/PBSCM_PGF/PB_SCM_PGF.py +++ b/cdmir/discovery/constraint/PBSCM_PGF/PB_SCM_PGF.py @@ -1,7 +1,7 @@ -from cdmir.effect.PBSCM_PGF.util import * +from cdmir.discovery.constraint.PBSCM_PGF.util import * from tqdm import tqdm from copy import deepcopy -from cdmir.effect.PBSCM_PGF.CCARankTest import CCARankTester +from cdmir.discovery.constraint.PBSCM_PGF.CCARankTest import CCARankTester np.set_printoptions(linewidth=np.inf, precision=3) diff --git a/cdmir/discovery/SHP/__init__.py b/cdmir/discovery/constraint/PBSCM_PGF/__init__.py similarity index 100% rename from cdmir/discovery/SHP/__init__.py rename to cdmir/discovery/constraint/PBSCM_PGF/__init__.py diff --git a/cdmir/effect/PBSCM_PGF/util.py b/cdmir/discovery/constraint/PBSCM_PGF/util.py similarity index 100% rename from cdmir/effect/PBSCM_PGF/util.py rename to cdmir/discovery/constraint/PBSCM_PGF/util.py diff --git a/cdmir/discovery/SHP/Generate_Hawkes_data_from_tick.py b/cdmir/discovery/funtional_based/SHP/Generate_Hawkes_data_from_tick.py similarity index 97% rename from cdmir/discovery/SHP/Generate_Hawkes_data_from_tick.py rename to cdmir/discovery/funtional_based/SHP/Generate_Hawkes_data_from_tick.py index 9c30465..7aa1260 100644 --- a/cdmir/discovery/SHP/Generate_Hawkes_data_from_tick.py +++ b/cdmir/discovery/funtional_based/SHP/Generate_Hawkes_data_from_tick.py @@ -1,109 +1,109 @@ -import pandas as pd -from itertools import product -import numpy as np -from tick.hawkes import SimuHawkes, HawkesKernelExp -import argparse -import os - -print(os.getcwd()) - - -def check_and_create(path: str): - if os.path.exists(path): - return - os.mkdir(path) - - -parser = argparse.ArgumentParser(description="data") -parser.add_argument('-sd', '--save_dir', type=str, help='save dir', default='./data/Hawkes_data_from_tick') -parser.add_argument('-et', '--exp_tag', type=int, help='data type', default=0) - - -def get_artificial_data(mu_range, alpha_range, n, sample_size=30000, out_degree_rate=1.5, NE_num=40, decay=0.1, - time_interval=None, seed=None): - print( - f'n={n},mu_range={mu_range},alpha_range={alpha_range},edge_num={round(out_degree_rate * n)},sample_size={sample_size}') - rand_state = np.random.RandomState(seed=seed) - - edge_mat = np.zeros([n, n]) - edge_select = list(filter(lambda i: i[0] < i[1], product(range(n), range(n)))) - rand_state.shuffle(edge_select) - for edge_ind in edge_select[:round(out_degree_rate * n)]: - edge_mat[edge_ind] = 1 - mu = rand_state.uniform(*mu_range, n) - - alpha = rand_state.uniform(*alpha_range, [n, n]) - alpha = edge_mat * alpha - - hawkes = SimuHawkes(baseline=mu, max_jumps=sample_size / NE_num, verbose=False, seed=seed) - for i in range(n): - for j in range(n): - if (alpha[i, j] == 0): - continue - hawkes.set_kernel(j, i, HawkesKernelExp(alpha[i, j], decay)) - - event_dict = dict() - for node in range(NE_num): - hawkes.reset() - hawkes.simulate() - event_dict[node] = hawkes.timestamps - - event_list = [] - - - for node in event_dict: - for event_name in range(n): - for timestamp in event_dict[node][event_name]: - event_list.append([node, timestamp, event_name]) - columns = ['seq_id', 'time_stamp', 'event_type'] - event_table = pd.DataFrame(event_list, columns=columns) - - if time_interval is not None: - event_table['time_stamp'] = (event_table['time_stamp'] / time_interval).astype('int') * time_interval - - events = [[event_table[(event_table['event_type'] == i) & (event_table['seq_id'] == j)][ - 'time_stamp'].values.astype('float') for i in - np.sort((event_table['event_type']).unique())] for j in (event_table['seq_id']).unique()] - return event_table, edge_mat, alpha, mu, events - - -def INSEM_data(sample_size=10000, lambda_x=1, theta=0.5, lambda_e=1, seed=None): - rand_state = np.random.RandomState(seed=seed) - X = rand_state.poisson(lambda_x, sample_size) - Y = np.zeros(sample_size, dtype='int') - - def operator(X): - return sum(rand_state.binomial(1, theta, X)) - - for i in range(sample_size): - Y[i] = operator(X[i]) - Y = Y + rand_state.poisson(lambda_e, sample_size) - t = 0 - df = pd.DataFrame() - for i, n in enumerate(X): - term = [] - for j in range(n): - term.append([0, t, 0]) - t = t + 1 - df = pd.concat([df, pd.DataFrame(term)]) - t = 0 - for i, n in enumerate(Y): - term = [] - for j in range(n): - term.append([0, t, 1]) - t = t + 1 - df = pd.concat([df, pd.DataFrame(term)]) - df.columns = ['seq_id', 'time_stamp', 'event_type'] - return df - - -def generate_data(n, mu_range_str, alpha_range_str, sample_size, out_degree_rate, NE_num, decay, seed=None): - alpha_range = tuple([float(i) for i in alpha_range_str.split(',')]) - mu_range = tuple([float(i) for i in mu_range_str.split(',')]) - event_table, edge_mat, alpha, mu, events = get_artificial_data(mu_range, alpha_range, n, sample_size=sample_size, - out_degree_rate=out_degree_rate, NE_num=NE_num, - decay=decay, seed=seed) - - return event_table, edge_mat, alpha, mu, events - - +import pandas as pd +from itertools import product +import numpy as np +from tick.hawkes import SimuHawkes, HawkesKernelExp +import argparse +import os + +print(os.getcwd()) + + +def check_and_create(path: str): + if os.path.exists(path): + return + os.mkdir(path) + + +parser = argparse.ArgumentParser(description="data") +parser.add_argument('-sd', '--save_dir', type=str, help='save dir', default='./data/Hawkes_data_from_tick') +parser.add_argument('-et', '--exp_tag', type=int, help='data type', default=0) + + +def get_artificial_data(mu_range, alpha_range, n, sample_size=30000, out_degree_rate=1.5, NE_num=40, decay=0.1, + time_interval=None, seed=None): + print( + f'n={n},mu_range={mu_range},alpha_range={alpha_range},edge_num={round(out_degree_rate * n)},sample_size={sample_size}') + rand_state = np.random.RandomState(seed=seed) + + edge_mat = np.zeros([n, n]) + edge_select = list(filter(lambda i: i[0] < i[1], product(range(n), range(n)))) + rand_state.shuffle(edge_select) + for edge_ind in edge_select[:round(out_degree_rate * n)]: + edge_mat[edge_ind] = 1 + mu = rand_state.uniform(*mu_range, n) + + alpha = rand_state.uniform(*alpha_range, [n, n]) + alpha = edge_mat * alpha + + hawkes = SimuHawkes(baseline=mu, max_jumps=sample_size / NE_num, verbose=False, seed=seed) + for i in range(n): + for j in range(n): + if (alpha[i, j] == 0): + continue + hawkes.set_kernel(j, i, HawkesKernelExp(alpha[i, j], decay)) + + event_dict = dict() + for node in range(NE_num): + hawkes.reset() + hawkes.simulate() + event_dict[node] = hawkes.timestamps + + event_list = [] + + + for node in event_dict: + for event_name in range(n): + for timestamp in event_dict[node][event_name]: + event_list.append([node, timestamp, event_name]) + columns = ['seq_id', 'time_stamp', 'event_type'] + event_table = pd.DataFrame(event_list, columns=columns) + + if time_interval is not None: + event_table['time_stamp'] = (event_table['time_stamp'] / time_interval).astype('int') * time_interval + + events = [[event_table[(event_table['event_type'] == i) & (event_table['seq_id'] == j)][ + 'time_stamp'].values.astype('float') for i in + np.sort((event_table['event_type']).unique())] for j in (event_table['seq_id']).unique()] + return event_table, edge_mat, alpha, mu, events + + +def INSEM_data(sample_size=10000, lambda_x=1, theta=0.5, lambda_e=1, seed=None): + rand_state = np.random.RandomState(seed=seed) + X = rand_state.poisson(lambda_x, sample_size) + Y = np.zeros(sample_size, dtype='int') + + def operator(X): + return sum(rand_state.binomial(1, theta, X)) + + for i in range(sample_size): + Y[i] = operator(X[i]) + Y = Y + rand_state.poisson(lambda_e, sample_size) + t = 0 + df = pd.DataFrame() + for i, n in enumerate(X): + term = [] + for j in range(n): + term.append([0, t, 0]) + t = t + 1 + df = pd.concat([df, pd.DataFrame(term)]) + t = 0 + for i, n in enumerate(Y): + term = [] + for j in range(n): + term.append([0, t, 1]) + t = t + 1 + df = pd.concat([df, pd.DataFrame(term)]) + df.columns = ['seq_id', 'time_stamp', 'event_type'] + return df + + +def generate_data(n, mu_range_str, alpha_range_str, sample_size, out_degree_rate, NE_num, decay, seed=None): + alpha_range = tuple([float(i) for i in alpha_range_str.split(',')]) + mu_range = tuple([float(i) for i in mu_range_str.split(',')]) + event_table, edge_mat, alpha, mu, events = get_artificial_data(mu_range, alpha_range, n, sample_size=sample_size, + out_degree_rate=out_degree_rate, NE_num=NE_num, + decay=decay, seed=seed) + + return event_table, edge_mat, alpha, mu, events + + diff --git a/cdmir/discovery/SHP/SHP.py b/cdmir/discovery/funtional_based/SHP/SHP.py similarity index 97% rename from cdmir/discovery/SHP/SHP.py rename to cdmir/discovery/funtional_based/SHP/SHP.py index 4c70b9b..7eb7c09 100644 --- a/cdmir/discovery/SHP/SHP.py +++ b/cdmir/discovery/funtional_based/SHP/SHP.py @@ -1,455 +1,455 @@ -from __future__ import print_function -import pandas as pd -import numpy as np -from tqdm import tqdm -import networkx as nx -from itertools import count, product -from Generate_Hawkes_data_from_tick import generate_data -from functools import partialmethod - - -#tqdm.__init__ = partialmethod(tqdm.__init__, disable=True) - -__MIN__ = -np.inf - - -def get_interval_events(event_table, time_interval=20): - """ - Discretize event timestamps into time intervals and organize events by sequence ID and event type - - :param event_table : pandas.DataFrame - Event data table containing columns 'time_stamp', 'event_type', and 'seq_id' - :param time_interval : int, default=20 - Time interval size for discretizing timestamps - - - :return: list - 3D list where the first dimension is sequence ID, - second dimension is event type, - and third dimension is array of timestamps for that sequence and event type - """ - - # Discretize timestamps to specified time intervals - event_table['time_stamp'] = (event_table['time_stamp'] / time_interval).astype('int') * time_interval - - # Group and organize timestamps by sequence ID and event type - events = [[event_table[(event_table['event_type'] == i) & (event_table['seq_id'] == j)][ - 'time_stamp'].values.astype('float') for i in - np.sort((event_table['event_type']).unique())] for j in tqdm((event_table['seq_id']).unique())] - return events - - -def check_DAG(edge_mat): - c_g = nx.from_numpy_array(edge_mat - np.diag(np.diag(edge_mat)), create_using=nx.DiGraph) - return nx.is_directed_acyclic_graph(c_g) - - -class SHP(object): - """ - Python implementation of SHP algorithm - - Reference: - ---------- - [1] Qiao J, Cai R, Wu S, et al. Structural hawkes processes for learning causal structure from discrete-time event sequences[J]. arXiv preprint arXiv:2305.05986, 2023. - """ - def __init__(self, event_table: pd.DataFrame, decay, time_interval=None, - init_structure: np.array = None, - penalty='BIC', seed=None, reg=3.0): - """ - Construct the SHP model - :param event_table: A pandas.DataFrame of events with columns ['seq_id', 'time_stamp', 'event_type'] - :param decay: The decay used in the exponential kernel - :param init_structure: adj of causal structure of prior knowledge - :param penalty: 'BIC' or 'AIC' penalty - """ - - self.random_state = np.random.RandomState(seed) - self.reg = reg - self.time_interval = time_interval # Delta t - - # create the event_table - self.event_table, self.event_names = self.get_event_table(event_table) - - # set the start timestamp to zero - for seq_id in np.unique(self.event_table["seq_id"].values): - seq_index = self.event_table[self.event_table["seq_id"] == seq_id].index.tolist() - self.event_table.loc[seq_index, "max_time_stamp"] = \ - self.event_table.loc[seq_index]["max_time_stamp"] - self.event_table.loc[seq_index]["time_stamp"].min() - self.event_table.loc[seq_index, "time_stamp"] = \ - self.event_table.loc[seq_index]["time_stamp"] - self.event_table.loc[seq_index]["time_stamp"].min() - - # store the calculated likelihood - self.hist_likelihood = dict() - for i in range(len(self.event_names)): - self.hist_likelihood[i] = dict() - - if penalty not in {'BIC', 'AIC'}: - raise Exception('Penalty is not supported') - self.penalty = penalty - - self.decay = decay # the decay coefficient of kappa function - self.n = len(self.event_names) # num of event type - self.T = self.event_table.groupby('seq_id').apply(lambda i: (i['time_stamp'].max())).sum() # total time span - self.T_each_seq = self.event_table.groupby('seq_id'). \ - apply(lambda i: (i['time_stamp'].max())) # the last moment of each event sequence - - # Initializing Structs - if init_structure is None: - self.init_structure = np.zeros([self.n, self.n]) - elif not ((init_structure == 1) | (init_structure == 0)).all(): - raise ValueError('Elements of the adjacency matrix need to be 0 or 1') - else: - self.init_structure = np.array(init_structure) - - X_dict = dict() - for seq_id, time_stamp, _, times, type_ind, _ in self.event_table.values: - if (seq_id, time_stamp) not in X_dict: - X_dict[(seq_id, time_stamp)] = [0] * self.n - X_dict[(seq_id, time_stamp)][type_ind] = times - self.X_df = pd.DataFrame(X_dict).T - - self.X = self.X_df.values - self.sum_t_X_kappa, self.decay_effect_integral_to_T = self.calculate_influence_of_each_event() - - def calculate_influence_of_each_event(self): - """ - calculate the influence of each event - - :return: - """ - sum_t_X_kappa = np.zeros_like(self.X, dtype='float64') - decay_effect_integral_to_T = self.X_df.copy() - - for ind, (seq_id, time_stamp) in tqdm(enumerate(self.X_df.index)): - # calculate the integral of decay function on time - decay_effect_integral_to_T.iloc[ind] = \ - self.X_df.iloc[ind] * ((1 - np.exp(-self.decay * (self.T_each_seq[seq_id] - time_stamp))) / self.decay) - - start_ind = ind - start_seq_id, start_time_stamp = self.X_df.index[start_ind] - # the influence on subsequent timestamp when the event occurs - next_ind = start_ind - while start_seq_id == seq_id: # the influence only spread on the same sequence - kap = self.kappa(start_time_stamp - time_stamp) - if kap < 0.0001: - break - - X_kappa = self.X[ind] * kap - sum_t_X_kappa[next_ind] += X_kappa # record the influence - next_ind += 1 - if next_ind >= len(self.X): - break - start_seq_id, start_time_stamp = self.X_df.index[next_ind] - return sum_t_X_kappa, decay_effect_integral_to_T - - # decay function - def kappa(self, t): - y = np.exp(-self.decay * t) - return y - - # transfer event table from continuous time domain to the discrete time domain - def get_event_table(self, event_table: pd.DataFrame): - event_table = event_table.copy() - event_table.columns = ['seq_id', 'time_stamp', 'event_type'] - if self.time_interval is not None: - event_table['time_stamp'] = (event_table['time_stamp'] / self.time_interval).astype( - 'int') * self.time_interval - - event_table['times'] = np.zeros(len(event_table)) - event_table = event_table.groupby(['seq_id', 'time_stamp', 'event_type']).count().reset_index() - - event_ind = event_table['event_type'].astype('category') - event_table['type_ind'] = event_ind.cat.codes - event_names = event_ind.cat.categories - - max_time = event_table.groupby('seq_id').apply(lambda i: i['time_stamp'].max()) - event_table = pd.merge(event_table, pd.DataFrame(max_time, columns=['max_time_stamp']).reset_index()) - - event_table.sort_values(['seq_id', 'time_stamp', 'type_ind']) - - return event_table, event_names - - def EM(self, edge_mat): - """ - EM (Expectation-Maximization) module with using hill climb - :param edge_mat: Adjacency matrix - :return: (likelihood, alpha matrix, mu vector) - """ - - if not check_DAG(edge_mat): - return __MIN__, np.zeros([len(self.event_names), len(self.event_names)]), np.zeros( - len(self.event_names)) - - # Initialize alpha and mu parameters - alpha = self.random_state.uniform(0, 1, [len(self.event_names), len(self.event_names)]) - alpha = alpha * edge_mat - mu = np.ones(len(self.event_names)) - L = 0 - - # calculate the likelihood for each event type i - for i in (range(len(self.event_names))): - # Get the parent node set of event i (i.e., all events pointing to i - Pa_i = tuple(np.where(edge_mat[:, i] == 1)[0]) - - try: - # Try to retrieve previously calculated results from cache - Li = self.hist_likelihood[i][tuple(Pa_i)][0] - mu[i] = self.hist_likelihood[i][tuple(Pa_i)][2] - for j in Pa_i: - alpha[j, i] = self.hist_likelihood[i][tuple(Pa_i)][1][j] - L += Li - except Exception as e: - # If not cached, optimize parameters iteratively using EM algorithm - Li = __MIN__ - - while 1: - # the first term of likelihood function - lambda_for_i = (self.sum_t_X_kappa * alpha[:, i]).sum(1) + mu[i] - # the second term of likelihood function - X_log_lambda = (self.X[:, i] * np.log(lambda_for_i)).sum() - lambda_i_sum = (((1 / self.decay) * self.X).sum(0) * alpha[:, i].T).sum() + mu[i] * self.T - - # calculate the likelihood - new_Li = -lambda_i_sum + X_log_lambda - - # Iteration termination condition - gain = new_Li - Li - if gain < 0.0085: - Li = new_Li - L += Li - Pa_i_alpha = dict() - for j in Pa_i: - Pa_i_alpha[j] = alpha[j, i] - self.hist_likelihood[i][tuple(Pa_i)] = (Li, Pa_i_alpha, mu[i]) - break - Li = new_Li - - # update mu - mu[i] = ((mu[i] / lambda_for_i) * self.X[:, i]).sum() / self.T - # update alpha - for j in Pa_i: - q_alpha = alpha[j, i] * self.sum_t_X_kappa[:, j] / lambda_for_i - upper = (q_alpha * self.X[:, i]).sum() - - lower = self.decay_effect_integral_to_T.sum(0)[j] * self.time_interval - if lower == 0: - alpha[j, i] = 0 - continue - alpha[j, i] = upper / lower - - i += 1 - # Adjust likelihood value according to specified penalty method (AIC or BIC - if self.penalty == 'AIC': - return L - (len(self.event_names) + edge_mat.sum()), alpha, mu - if self.penalty == 'BIC': - return L - (len(self.event_names) + edge_mat.sum()) * np.log( - self.event_table['times'].sum()) * self.reg, alpha, mu - - - def EM_not_HC(self, edge_mat): - """ - EM (Expectation-Maximization) module without using hill climb - :param edge_mat: Adjacency matrix - - :return: (likelihood, alpha matrix, mu vector) - """ - - # Initialize alpha and mu parameters - alpha = self.random_state.uniform(0, 1, [len(self.event_names), len(self.event_names)]) - alpha = alpha * edge_mat - mu = np.ones(len(self.event_names)) - L = 0 - - # Perform parameter estimation for each event i - for i in (range(len(self.event_names))): - # Get the parent node set of event i (i.e., all events pointing to i) - Pa_i = set(np.where(edge_mat[:, i] == 1)[0]) - - try: - # Try to retrieve previously calculated results from cache - Li = self.hist_likelihood[i][tuple(Pa_i)][0] - mu[i] = self.hist_likelihood[i][tuple(Pa_i)][2] - for j in Pa_i: - alpha[j, i] = self.hist_likelihood[i][tuple(Pa_i)][1][j] - L += Li - except Exception as e: - # If not cached, optimize parameters iteratively using EM algorithm - Li = __MIN__ - - while 1: - # the first term of likelihood function - lambda_for_i = (self.sum_t_X_kappa * alpha[:, i]).sum(1) + mu[i] - # the second term of likelihood function - X_log_lambda = (self.X[:, i] * np.log(lambda_for_i)).sum() - lambda_i_sum = (((1 / self.decay) * self.X).sum(0) * alpha[:, i].T).sum() + mu[i] * self.T - - # Current likelihood value - new_Li = -lambda_i_sum + X_log_lambda - # Iteration termination condition - gain = new_Li - Li - - # Check if termination condition is met - if gain <= 0.0085: - Li = new_Li - L += Li - # Cache current parameter results - Pa_i_alpha = dict() - for j in Pa_i: - Pa_i_alpha[j] = alpha[j, i] - self.hist_likelihood[i][tuple(Pa_i)] = (Li, Pa_i_alpha, mu[i]) - break - Li = new_Li - - # update mu - mu[i] = ((mu[i] / lambda_for_i) * self.X[:, i]).sum() / self.T - # update alpha - for j in Pa_i: - q_alpha = alpha[j, i] * self.sum_t_X_kappa[:, j] / lambda_for_i - upper = (q_alpha * self.X[:, i]).sum() - - lower = self.decay_effect_integral_to_T.sum(0)[j] * self.time_interval - if lower == 0: - alpha[j, i] = 0 - continue - alpha[j, i] = upper / lower - - i += 1 - - # Adjust likelihood value according to specified penalty method (AIC or BIC) - if self.penalty == 'AIC': - return L - (len(self.event_names) + edge_mat.sum())* self.reg, alpha, mu - if self.penalty == 'BIC': - return L - (len(self.event_names) + edge_mat.sum()) * np.log( - self.event_table['times'].sum()) * self.reg, alpha, mu - - # the searching module for new edges - def one_step_change_iterator(self, edge_mat): - return map(lambda e: self.one_step_change(edge_mat, e), - product(range(len(self.event_names)), - range(len(self.event_names)), range(3))) - - def one_step_change(self, edge_mat, e): - j, i = e[0], e[1] - if j == i: - return edge_mat - new_graph = edge_mat.copy() - if e[2] == 0: - new_graph[j, i] = 0 - new_graph[i, j] = 0 - return new_graph - elif e[2] == 1: - new_graph[j, i] = 1 - new_graph[i, j] = 0 - return new_graph - else: - new_graph[j, i] = 0 - new_graph[i, j] = 1 - return new_graph - - def Hill_Climb(self): - - # Initialize the adjacency matrix - edge_mat = self.init_structure - result = self.EM(edge_mat) - - L = result[0] - while 1: - stop_tag = True - for new_edge_mat in tqdm(list(self.one_step_change_iterator(edge_mat)), mininterval=5): - new_result = self.EM(new_edge_mat) - new_L = new_result[0] - # Termination condition: no adjacency matrix with higher likelihood appears - if new_L > L: - result = new_result - L = new_L - # if there is a new edge can be added, then set the stop_tag=False and continue searching - stop_tag = False - edge_mat = new_edge_mat - - if stop_tag: - return result - - def train_model(self,hill_climb=True): - """ - Train the model and return the result - :param hill_climb: use hill_climb or not - :return:(likelihood,alpha,mu) - """ - if hill_climb: - res = self.Hill_Climb() - else: - res = self.EM_not_HC(np.ones([self.n, self.n]) - np.eye(self.n, self.n)) - return res - - -def SHP_exp(n, sample_size, out_degree_rate, mu_range_str, alpha_range_str, decay, penalty='BIC', - NE_num=40, model_decay=0.35, seed=0, - time_interval=5, hill_climb=True, reg=0.85): - - """ - Execute SHP model experiment, generate simulated data and perform parameter estimation - - :param n: Number of nodes - :param sample_size: Sample size - :param out_degree_rate: Average out-degree ratio of nodes - :param mu_range_str: String representing the range of mu parameter values - :param alpha_range_str: String representing the range of alpha parameter values - :param decay: Decay parameter - :param penalty: Type of penalty term, default is 'BIC' - :param NE_num: Number of non-events, default is 40 - :param model_decay: Model decay parameter, default is 0.35 - :param seed: Random seed, default is 0 - :param time_interval: Time interval, default is 5 - :param hill_climb: Whether to use hill climbing algorithm, default is True - :param reg: Regularization parameter, default is 0.85 - - :return: Tuple containing: - - likelihood (float): Likelihood value - - fited_alpha (array): Fitted alpha parameters - - fited_mu (array): Fitted mu parameters - - real_edge_mat (array): Real adjacency matrix - - real_alpha (array): Real alpha parameters - - real_mu (array): Real mu parameters - """ - - - # Generate simulated data, including event table, real adjacency matrix and parameters - event_table, real_edge_mat, real_alpha, real_mu, events = generate_data(n=n, sample_size=sample_size, - out_degree_rate=out_degree_rate, - mu_range_str=mu_range_str, - alpha_range_str=alpha_range_str, - NE_num=NE_num, decay=decay, - seed=seed) - - # Build model parameter dictionary - param_dict = { - "decay": model_decay, - "reg": reg, - "time_interval": time_interval, - "penalty": penalty - } - # Initialize the SHP model instance - self = SHP(event_table, **param_dict) - - # Select optimization algorithm for parameter estimation based on parameters - res1=self.train() - - # Extract estimation results - likelihood = res1[0] - fited_alpha = res1[1] - fited_mu = res1[2] - return likelihood, fited_alpha, fited_mu, real_edge_mat, real_alpha, real_mu - - -if __name__=='__main__': - from utils import get_performance - - likelihood, fited_alpha, fited_mu, real_edge_mat, real_alpha, real_mu = SHP_exp(n=20, sample_size=20000, - out_degree_rate=1.5, - mu_range_str="0.00005,0.0001", - alpha_range_str="0.5,0.7", - decay=5, model_decay=0.35, seed=0, - time_interval=5, penalty='BIC', - hill_climb=False, reg=0.85) - res = get_performance(fited_alpha, real_edge_mat) - print(res) +from __future__ import print_function +import pandas as pd +import numpy as np +from tqdm import tqdm +import networkx as nx +from itertools import count, product +from Generate_Hawkes_data_from_tick import generate_data +from functools import partialmethod + + +#tqdm.__init__ = partialmethod(tqdm.__init__, disable=True) + +__MIN__ = -np.inf + + +def get_interval_events(event_table, time_interval=20): + """ + Discretize event timestamps into time intervals and organize events by sequence ID and event type + + :param event_table : pandas.DataFrame + Event data table containing columns 'time_stamp', 'event_type', and 'seq_id' + :param time_interval : int, default=20 + Time interval size for discretizing timestamps + + + :return: list + 3D list where the first dimension is sequence ID, + second dimension is event type, + and third dimension is array of timestamps for that sequence and event type + """ + + # Discretize timestamps to specified time intervals + event_table['time_stamp'] = (event_table['time_stamp'] / time_interval).astype('int') * time_interval + + # Group and organize timestamps by sequence ID and event type + events = [[event_table[(event_table['event_type'] == i) & (event_table['seq_id'] == j)][ + 'time_stamp'].values.astype('float') for i in + np.sort((event_table['event_type']).unique())] for j in tqdm((event_table['seq_id']).unique())] + return events + + +def check_DAG(edge_mat): + c_g = nx.from_numpy_array(edge_mat - np.diag(np.diag(edge_mat)), create_using=nx.DiGraph) + return nx.is_directed_acyclic_graph(c_g) + + +class SHP(object): + """ + Python implementation of SHP algorithm + + Reference: + ---------- + [1] Qiao J, Cai R, Wu S, et al. Structural hawkes processes for learning causal structure from discrete-time event sequences[J]. arXiv preprint arXiv:2305.05986, 2023. + """ + def __init__(self, event_table: pd.DataFrame, decay, time_interval=None, + init_structure: np.array = None, + penalty='BIC', seed=None, reg=3.0): + """ + Construct the SHP model + :param event_table: A pandas.DataFrame of events with columns ['seq_id', 'time_stamp', 'event_type'] + :param decay: The decay used in the exponential kernel + :param init_structure: adj of causal structure of prior knowledge + :param penalty: 'BIC' or 'AIC' penalty + """ + + self.random_state = np.random.RandomState(seed) + self.reg = reg + self.time_interval = time_interval # Delta t + + # create the event_table + self.event_table, self.event_names = self.get_event_table(event_table) + + # set the start timestamp to zero + for seq_id in np.unique(self.event_table["seq_id"].values): + seq_index = self.event_table[self.event_table["seq_id"] == seq_id].index.tolist() + self.event_table.loc[seq_index, "max_time_stamp"] = \ + self.event_table.loc[seq_index]["max_time_stamp"] - self.event_table.loc[seq_index]["time_stamp"].min() + self.event_table.loc[seq_index, "time_stamp"] = \ + self.event_table.loc[seq_index]["time_stamp"] - self.event_table.loc[seq_index]["time_stamp"].min() + + # store the calculated likelihood + self.hist_likelihood = dict() + for i in range(len(self.event_names)): + self.hist_likelihood[i] = dict() + + if penalty not in {'BIC', 'AIC'}: + raise Exception('Penalty is not supported') + self.penalty = penalty + + self.decay = decay # the decay coefficient of kappa function + self.n = len(self.event_names) # num of event type + self.T = self.event_table.groupby('seq_id').apply(lambda i: (i['time_stamp'].max())).sum() # total time span + self.T_each_seq = self.event_table.groupby('seq_id'). \ + apply(lambda i: (i['time_stamp'].max())) # the last moment of each event sequence + + # Initializing Structs + if init_structure is None: + self.init_structure = np.zeros([self.n, self.n]) + elif not ((init_structure == 1) | (init_structure == 0)).all(): + raise ValueError('Elements of the adjacency matrix need to be 0 or 1') + else: + self.init_structure = np.array(init_structure) + + X_dict = dict() + for seq_id, time_stamp, _, times, type_ind, _ in self.event_table.values: + if (seq_id, time_stamp) not in X_dict: + X_dict[(seq_id, time_stamp)] = [0] * self.n + X_dict[(seq_id, time_stamp)][type_ind] = times + self.X_df = pd.DataFrame(X_dict).T + + self.X = self.X_df.values + self.sum_t_X_kappa, self.decay_effect_integral_to_T = self.calculate_influence_of_each_event() + + def calculate_influence_of_each_event(self): + """ + calculate the influence of each event + + :return: + """ + sum_t_X_kappa = np.zeros_like(self.X, dtype='float64') + decay_effect_integral_to_T = self.X_df.copy() + + for ind, (seq_id, time_stamp) in tqdm(enumerate(self.X_df.index)): + # calculate the integral of decay function on time + decay_effect_integral_to_T.iloc[ind] = \ + self.X_df.iloc[ind] * ((1 - np.exp(-self.decay * (self.T_each_seq[seq_id] - time_stamp))) / self.decay) + + start_ind = ind + start_seq_id, start_time_stamp = self.X_df.index[start_ind] + # the influence on subsequent timestamp when the event occurs + next_ind = start_ind + while start_seq_id == seq_id: # the influence only spread on the same sequence + kap = self.kappa(start_time_stamp - time_stamp) + if kap < 0.0001: + break + + X_kappa = self.X[ind] * kap + sum_t_X_kappa[next_ind] += X_kappa # record the influence + next_ind += 1 + if next_ind >= len(self.X): + break + start_seq_id, start_time_stamp = self.X_df.index[next_ind] + return sum_t_X_kappa, decay_effect_integral_to_T + + # decay function + def kappa(self, t): + y = np.exp(-self.decay * t) + return y + + # transfer event table from continuous time domain to the discrete time domain + def get_event_table(self, event_table: pd.DataFrame): + event_table = event_table.copy() + event_table.columns = ['seq_id', 'time_stamp', 'event_type'] + if self.time_interval is not None: + event_table['time_stamp'] = (event_table['time_stamp'] / self.time_interval).astype( + 'int') * self.time_interval + + event_table['times'] = np.zeros(len(event_table)) + event_table = event_table.groupby(['seq_id', 'time_stamp', 'event_type']).count().reset_index() + + event_ind = event_table['event_type'].astype('category') + event_table['type_ind'] = event_ind.cat.codes + event_names = event_ind.cat.categories + + max_time = event_table.groupby('seq_id').apply(lambda i: i['time_stamp'].max()) + event_table = pd.merge(event_table, pd.DataFrame(max_time, columns=['max_time_stamp']).reset_index()) + + event_table.sort_values(['seq_id', 'time_stamp', 'type_ind']) + + return event_table, event_names + + def EM(self, edge_mat): + """ + EM (Expectation-Maximization) module with using hill climb + :param edge_mat: Adjacency matrix + :return: (likelihood, alpha matrix, mu vector) + """ + + if not check_DAG(edge_mat): + return __MIN__, np.zeros([len(self.event_names), len(self.event_names)]), np.zeros( + len(self.event_names)) + + # Initialize alpha and mu parameters + alpha = self.random_state.uniform(0, 1, [len(self.event_names), len(self.event_names)]) + alpha = alpha * edge_mat + mu = np.ones(len(self.event_names)) + L = 0 + + # calculate the likelihood for each event type i + for i in (range(len(self.event_names))): + # Get the parent node set of event i (i.e., all events pointing to i + Pa_i = tuple(np.where(edge_mat[:, i] == 1)[0]) + + try: + # Try to retrieve previously calculated results from cache + Li = self.hist_likelihood[i][tuple(Pa_i)][0] + mu[i] = self.hist_likelihood[i][tuple(Pa_i)][2] + for j in Pa_i: + alpha[j, i] = self.hist_likelihood[i][tuple(Pa_i)][1][j] + L += Li + except Exception as e: + # If not cached, optimize parameters iteratively using EM algorithm + Li = __MIN__ + + while 1: + # the first term of likelihood function + lambda_for_i = (self.sum_t_X_kappa * alpha[:, i]).sum(1) + mu[i] + # the second term of likelihood function + X_log_lambda = (self.X[:, i] * np.log(lambda_for_i)).sum() + lambda_i_sum = (((1 / self.decay) * self.X).sum(0) * alpha[:, i].T).sum() + mu[i] * self.T + + # calculate the likelihood + new_Li = -lambda_i_sum + X_log_lambda + + # Iteration termination condition + gain = new_Li - Li + if gain < 0.0085: + Li = new_Li + L += Li + Pa_i_alpha = dict() + for j in Pa_i: + Pa_i_alpha[j] = alpha[j, i] + self.hist_likelihood[i][tuple(Pa_i)] = (Li, Pa_i_alpha, mu[i]) + break + Li = new_Li + + # update mu + mu[i] = ((mu[i] / lambda_for_i) * self.X[:, i]).sum() / self.T + # update alpha + for j in Pa_i: + q_alpha = alpha[j, i] * self.sum_t_X_kappa[:, j] / lambda_for_i + upper = (q_alpha * self.X[:, i]).sum() + + lower = self.decay_effect_integral_to_T.sum(0)[j] * self.time_interval + if lower == 0: + alpha[j, i] = 0 + continue + alpha[j, i] = upper / lower + + i += 1 + # Adjust likelihood value according to specified penalty method (AIC or BIC + if self.penalty == 'AIC': + return L - (len(self.event_names) + edge_mat.sum()), alpha, mu + if self.penalty == 'BIC': + return L - (len(self.event_names) + edge_mat.sum()) * np.log( + self.event_table['times'].sum()) * self.reg, alpha, mu + + + def EM_not_HC(self, edge_mat): + """ + EM (Expectation-Maximization) module without using hill climb + :param edge_mat: Adjacency matrix + + :return: (likelihood, alpha matrix, mu vector) + """ + + # Initialize alpha and mu parameters + alpha = self.random_state.uniform(0, 1, [len(self.event_names), len(self.event_names)]) + alpha = alpha * edge_mat + mu = np.ones(len(self.event_names)) + L = 0 + + # Perform parameter estimation for each event i + for i in (range(len(self.event_names))): + # Get the parent node set of event i (i.e., all events pointing to i) + Pa_i = set(np.where(edge_mat[:, i] == 1)[0]) + + try: + # Try to retrieve previously calculated results from cache + Li = self.hist_likelihood[i][tuple(Pa_i)][0] + mu[i] = self.hist_likelihood[i][tuple(Pa_i)][2] + for j in Pa_i: + alpha[j, i] = self.hist_likelihood[i][tuple(Pa_i)][1][j] + L += Li + except Exception as e: + # If not cached, optimize parameters iteratively using EM algorithm + Li = __MIN__ + + while 1: + # the first term of likelihood function + lambda_for_i = (self.sum_t_X_kappa * alpha[:, i]).sum(1) + mu[i] + # the second term of likelihood function + X_log_lambda = (self.X[:, i] * np.log(lambda_for_i)).sum() + lambda_i_sum = (((1 / self.decay) * self.X).sum(0) * alpha[:, i].T).sum() + mu[i] * self.T + + # Current likelihood value + new_Li = -lambda_i_sum + X_log_lambda + # Iteration termination condition + gain = new_Li - Li + + # Check if termination condition is met + if gain <= 0.0085: + Li = new_Li + L += Li + # Cache current parameter results + Pa_i_alpha = dict() + for j in Pa_i: + Pa_i_alpha[j] = alpha[j, i] + self.hist_likelihood[i][tuple(Pa_i)] = (Li, Pa_i_alpha, mu[i]) + break + Li = new_Li + + # update mu + mu[i] = ((mu[i] / lambda_for_i) * self.X[:, i]).sum() / self.T + # update alpha + for j in Pa_i: + q_alpha = alpha[j, i] * self.sum_t_X_kappa[:, j] / lambda_for_i + upper = (q_alpha * self.X[:, i]).sum() + + lower = self.decay_effect_integral_to_T.sum(0)[j] * self.time_interval + if lower == 0: + alpha[j, i] = 0 + continue + alpha[j, i] = upper / lower + + i += 1 + + # Adjust likelihood value according to specified penalty method (AIC or BIC) + if self.penalty == 'AIC': + return L - (len(self.event_names) + edge_mat.sum())* self.reg, alpha, mu + if self.penalty == 'BIC': + return L - (len(self.event_names) + edge_mat.sum()) * np.log( + self.event_table['times'].sum()) * self.reg, alpha, mu + + # the searching module for new edges + def one_step_change_iterator(self, edge_mat): + return map(lambda e: self.one_step_change(edge_mat, e), + product(range(len(self.event_names)), + range(len(self.event_names)), range(3))) + + def one_step_change(self, edge_mat, e): + j, i = e[0], e[1] + if j == i: + return edge_mat + new_graph = edge_mat.copy() + if e[2] == 0: + new_graph[j, i] = 0 + new_graph[i, j] = 0 + return new_graph + elif e[2] == 1: + new_graph[j, i] = 1 + new_graph[i, j] = 0 + return new_graph + else: + new_graph[j, i] = 0 + new_graph[i, j] = 1 + return new_graph + + def Hill_Climb(self): + + # Initialize the adjacency matrix + edge_mat = self.init_structure + result = self.EM(edge_mat) + + L = result[0] + while 1: + stop_tag = True + for new_edge_mat in tqdm(list(self.one_step_change_iterator(edge_mat)), mininterval=5): + new_result = self.EM(new_edge_mat) + new_L = new_result[0] + # Termination condition: no adjacency matrix with higher likelihood appears + if new_L > L: + result = new_result + L = new_L + # if there is a new edge can be added, then set the stop_tag=False and continue searching + stop_tag = False + edge_mat = new_edge_mat + + if stop_tag: + return result + + def train_model(self,hill_climb=True): + """ + Train the model and return the result + :param hill_climb: use hill_climb or not + :return:(likelihood,alpha,mu) + """ + if hill_climb: + res = self.Hill_Climb() + else: + res = self.EM_not_HC(np.ones([self.n, self.n]) - np.eye(self.n, self.n)) + return res + + +def SHP_exp(n, sample_size, out_degree_rate, mu_range_str, alpha_range_str, decay, penalty='BIC', + NE_num=40, model_decay=0.35, seed=0, + time_interval=5, hill_climb=True, reg=0.85): + + """ + Execute SHP model experiment, generate simulated data and perform parameter estimation + + :param n: Number of nodes + :param sample_size: Sample size + :param out_degree_rate: Average out-degree ratio of nodes + :param mu_range_str: String representing the range of mu parameter values + :param alpha_range_str: String representing the range of alpha parameter values + :param decay: Decay parameter + :param penalty: Type of penalty term, default is 'BIC' + :param NE_num: Number of non-events, default is 40 + :param model_decay: Model decay parameter, default is 0.35 + :param seed: Random seed, default is 0 + :param time_interval: Time interval, default is 5 + :param hill_climb: Whether to use hill climbing algorithm, default is True + :param reg: Regularization parameter, default is 0.85 + + :return: Tuple containing: + - likelihood (float): Likelihood value + - fited_alpha (array): Fitted alpha parameters + - fited_mu (array): Fitted mu parameters + - real_edge_mat (array): Real adjacency matrix + - real_alpha (array): Real alpha parameters + - real_mu (array): Real mu parameters + """ + + + # Generate simulated data, including event table, real adjacency matrix and parameters + event_table, real_edge_mat, real_alpha, real_mu, events = generate_data(n=n, sample_size=sample_size, + out_degree_rate=out_degree_rate, + mu_range_str=mu_range_str, + alpha_range_str=alpha_range_str, + NE_num=NE_num, decay=decay, + seed=seed) + + # Build model parameter dictionary + param_dict = { + "decay": model_decay, + "reg": reg, + "time_interval": time_interval, + "penalty": penalty + } + # Initialize the SHP model instance + self = SHP(event_table, **param_dict) + + # Select optimization algorithm for parameter estimation based on parameters + res1=self.train() + + # Extract estimation results + likelihood = res1[0] + fited_alpha = res1[1] + fited_mu = res1[2] + return likelihood, fited_alpha, fited_mu, real_edge_mat, real_alpha, real_mu + + +if __name__=='__main__': + from utils import get_performance + + likelihood, fited_alpha, fited_mu, real_edge_mat, real_alpha, real_mu = SHP_exp(n=20, sample_size=20000, + out_degree_rate=1.5, + mu_range_str="0.00005,0.0001", + alpha_range_str="0.5,0.7", + decay=5, model_decay=0.35, seed=0, + time_interval=5, penalty='BIC', + hill_climb=False, reg=0.85) + res = get_performance(fited_alpha, real_edge_mat) + print(res) diff --git a/cdmir/effect/PBSCM_PGF/__init__.py b/cdmir/discovery/funtional_based/SHP/__init__.py similarity index 100% rename from cdmir/effect/PBSCM_PGF/__init__.py rename to cdmir/discovery/funtional_based/SHP/__init__.py diff --git a/docs/source/discovery_methods/constraint/PBSCM/pbscm.rst b/docs/source/discovery_methods/constraint/PBSCM/pbscm.rst new file mode 100644 index 0000000..e52ebcc --- /dev/null +++ b/docs/source/discovery_methods/constraint/PBSCM/pbscm.rst @@ -0,0 +1,47 @@ +PBSCM (Poisson Branching Structural Causal Model) +============================================= + +Introduction +------------ + +PBSCM is a functional-based causal discovery method that identifies causal relationships using high-order cumulants with path analysis. It assumes that the causal relationships can be modeled as a Poisson branching process. + + +Usage +----- + +.. code-block:: python + + + from cdmir.discovery.funtional_based.PBSCM.PB_SCM import PB_SCM + + # Initialize PBSCM model + pbscm = PB_SCM(data) + + # Get causal graph + causal_graph = pbscm.get_causal_graph(alpha=0.04, max_order=4, threshold=0) + + +Parameters +---------- + +**PB_SCM Class Initialization** + +- **data**: Input data matrix (n_samples x n_variables) + +**get_causal_graph** + + +- **alpha**: Confidence level (default: 0.04) +- **max_order**: The maximum order of Lambda_k (default: 4) +- **threshold**: Threshold when bootstrap test fails (default: 0) + +Returns +------- + +- **causal_graph**: Causal matrix of the data + +References +---------- + +.. [1] Qiao J, Xiang Y, Chen Z, et al. Causal discovery from poisson branching structural causal model using high-order cumulant with path analysis[C]//Proceedings of the AAAI Conference on Artificial Intelligence. 2024, 38(18): 20524-20531. diff --git a/docs/source/discovery_methods/constraint/PBSCM_PGF/pbscm_pgf.rst b/docs/source/discovery_methods/constraint/PBSCM_PGF/pbscm_pgf.rst new file mode 100644 index 0000000..8c38518 --- /dev/null +++ b/docs/source/discovery_methods/constraint/PBSCM_PGF/pbscm_pgf.rst @@ -0,0 +1,43 @@ +PBSCM_PGF (Poisson Branching Structural Causal Model using Probability Generating Function) +========================================================================================= + +Introduction +------------ + +PBSCM_PGF is a causal discovery method that identifies causal relationships using Probability Generating Functions (PGF). It extends the Poisson Branching Structural Causal Model (PBSCM) by leveraging PGFs to analyze the structure of causal relationships in data. + + +Usage +----- + +.. code-block:: python + + from cdmir.effect.PBSCM_PGF.PB_SCM_PGF import PBSCM_PGF + + # Initialize PBSCM_PGF model + pbscm_pgf = PBSCM_PGF(data) + + # Learn the causal graph + causal_graph = pbscm_pgf.get_causal_graph() + + + + +Parameters +---------- + +**PBSCM_PGF Class Initialization** + +- **data**: Input data matrix with shape (n_samples, n_features) + + +Returns +------- + +- **causal_graph**: Causal matrix of the data, where causal_graph[i][j] = 1 indicates a causal relationship from variable i to variable j, and 0 otherwise + + +References +---------- + +.. [1] Xiang Y, Qiao J, Liang Z, et al. On the identifiability of poisson branching structural causal model using probability generating function[J]. Advances in Neural Information Processing Systems, 2024, 37: 11664-11699. diff --git a/docs/source/discovery_methods/constraint/index.rst b/docs/source/discovery_methods/constraint/index.rst index 40c9798..e43d7c7 100644 --- a/docs/source/discovery_methods/constraint/index.rst +++ b/docs/source/discovery_methods/constraint/index.rst @@ -1,2 +1,11 @@ Constraint-based -================= \ No newline at end of file +================= +.. toctree:: + :maxdepth: 2 + + PBSCM/pbscm + +.. toctree:: + :maxdepth: 2 + + PBSCM_PGF/pbscm_pgf diff --git a/docs/source/discovery_methods/functional_based/ANM/anm.rst b/docs/source/discovery_methods/functional_based/ANM/anm.rst new file mode 100644 index 0000000..e7973b6 --- /dev/null +++ b/docs/source/discovery_methods/functional_based/ANM/anm.rst @@ -0,0 +1,37 @@ +ANM (Additive Noise Model) +============================ + +Introduction +------------ + +ANM is a functional-based causal discovery method that identifies causal relationships by testing the independence between the cause and the noise term in an additive noise model. + +Usage +----- + +.. code-block:: python + + from cdmir.discovery.funtional_based.anm.ANM import ANM + # Initialize ANM model + anm = ANM() + # Test causal direction + nonindepscore_forward, nonindepscore_backward = anm.cause_or_effect(X, Y) + + +Parameters +---------- + +- **x**: Input data array of shape (n,) or (n, 1) +- **y**: Output data array of shape (n,) or (n, 1) + +Returns +------- + +- **nonindepscore_forward**: HSIC statistic in the X→Y direction +- **nonindepscore_backward**: HSIC statistic in the Y→X direction + + +References +---------- + +.. [1] Hoyer, Patrik O., et al. "Nonlinear causal discovery with additive noise models." Advances in Neural Information Processing Systems. 2008. diff --git a/docs/source/discovery_methods/functional_based/SHP/shp.rst b/docs/source/discovery_methods/functional_based/SHP/shp.rst new file mode 100644 index 0000000..4de2b5c --- /dev/null +++ b/docs/source/discovery_methods/functional_based/SHP/shp.rst @@ -0,0 +1,56 @@ +SHP (Structural Hawkes Processes) +================================= + +Introduction +------------ + +SHP is a functional-based causal discovery method for learning causal structure from discrete-time event sequences. It extends the traditional Hawkes process by incorporating structural learning to identify causal relationships between event types. + +Usage +----- + +.. code-block:: python + + from cdmir.discovery.funtional_based.SHP.SHP import SHP + + # Initialize SHP model + shp = SHP( + event_table=event_table, + decay=0.35, + time_interval=5, + penalty='BIC', + reg=0.85, + seed=2025 + ) + + # Train the model with Hill Climb search + likelihood, alpha, mu = shp.train_model(hill_climb=True) + +Parameters +---------- + +**SHP Class Initialization** + +- **event_table**: Event data table containing columns 'seq_id', 'time_stamp', and 'event_type' +- **decay**: The decay coefficient of the exponential kernel +- **time_interval**: Time interval size for discretizing timestamps (default: None) +- **init_structure**: Initial adjacency matrix for causal structure (default: None) +- **penalty**: Penalty method, either 'BIC' or 'AIC' (default: 'BIC') +- **seed**: Random seed for reproducibility (default: None) +- **reg**: Regularization parameter (default: 3.0) + +**train_model Parameters** + +- **hill_climb**: Whether to use Hill Climb search for structure learning (default: True) + +Returns +------- + +- **likelihood**: Likelihood value of the learned model +- **alpha**: Matrix of causal effect parameters between event types, where :math:`\alpha_{ij}` represents the causal effect from event type :math:`i` to :math:`j` +- **mu**: Vector of base intensity parameters for each event type, where :math:`\mu_j` is the base intensity for event type :math:`j` + +References +---------- + +.. [1] Qiao J, Cai R, Wu S, et al. Structural hawkes processes for learning causal structure from discrete-time event sequences[J]. arXiv preprint arXiv:2305.05986, 2023. diff --git a/docs/source/discovery_methods/functional_based/index.rst b/docs/source/discovery_methods/functional_based/index.rst index b36ee41..23890c9 100644 --- a/docs/source/discovery_methods/functional_based/index.rst +++ b/docs/source/discovery_methods/functional_based/index.rst @@ -1,2 +1,12 @@ Functional-based -================= \ No newline at end of file +================= + +.. toctree:: + :maxdepth: 2 + + ANM/anm + +.. toctree:: + :maxdepth: 2 + + SHP/shp diff --git a/docs/source/effect_methods/index.rst b/docs/source/effect_methods/index.rst new file mode 100644 index 0000000..e15fd04 --- /dev/null +++ b/docs/source/effect_methods/index.rst @@ -0,0 +1,5 @@ +Discovery methods +================= + +In this section, we introduce effect methods implemented in CDMIR. + diff --git a/docs/source/index.rst b/docs/source/index.rst index 5c4627b..97a1dd5 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -27,6 +27,12 @@ Contents discovery_methods/index +.. toctree:: + :maxdepth: 1 + :caption: effect methods + + effect_methods/index + .. toctree:: :maxdepth: 1 :caption: Independence tests