diff --git a/main.py b/main.py index 64f0e54..71d6b20 100644 --- a/main.py +++ b/main.py @@ -1,324 +1,235 @@ import numpy as np import matplotlib.pyplot as plt import argparse import sys import pandas as pd from sklearn.metrics import mutual_info_score def soft_moduarity(wc, dc_inv, xc, lc): # Calculate the soft modularity as defined by equation (7) temp = np.linalg.multi_dot([dc_inv, xc, lc]) temp_t = np.transpose(temp) n = wc.shape[0] qc_s = np.trace(np.linalg.multi_dot([temp_t, wc, temp])) - np.linalg.multi_dot([np.ones((1, n)), temp, temp_t, wc, np.ones((n, 1))]) return qc_s def similarity(adj_mat, weighted, gamma=0.2): # calculate a similarity matrix from the adjacency matrix, basically an exponential scaling n, m = adj_mat.shape sim = np.zeros(adj_mat.shape) if weighted: for i in range(n): for j in range(m): if adj_mat[i, j] != 0: sim[i, j] = np.exp(-1.0/(gamma * adj_mat[i, j])) else: sim[i, j] = 0.0 else: sim = adj_mat / 2.0 for i in range(n): # populate the diagonal sim[i, i] = 1.0 sim /= np.sum(sim) return sim def normalize_rows(w: np.ndarray): # normalize matrix w row-wise - return w / np.linalg.norm(w, ord=1, axis=1, keepdims=True) def normalize_cols(xc: np.ndarray): # normalize matrix xc column-wise - return xc / np.linalg.norm(xc, ord=1, axis=0, keepdims=True) def kl_divergence(a_mat, b_mat): # Kullback-Leibler divergence for all non zero elements of a_mat and b_mat assert a_mat.shape == b_mat.shape n, m = a_mat.shape res = 0.0 for i in range(n): for j in range(m): if a_mat[i, j] > 1e-20 and b_mat[i, j] > 1e-20: res += a_mat[i, j] * np.log(a_mat[i, j] / b_mat[i, j]) - a_mat[i, j] + b_mat[i, j] return res -def param_update(yc, wc, alpha): +def param_update(yc, wc, alpha, is_plot=True): # Update function for xc (capital X and capital L), implementation of equations (4) and (5) n, m = yc.shape # xc_old = np.random.rand(xc.shape[0], xc.shape[1]) xc_old = np.random.rand(n, m) xc_old = normalize_cols(xc_old) # lc_old = lc_new lc_old = np.diag(np.random.rand(m)) lc_old /= np.trace(lc_old) - eps = 1e-6 + eps = 1e-5 max_iter = int(1e3) gamma = np.zeros(max_iter) for it in range(max_iter): wc_approx = xc_old.dot(lc_old.dot(xc_old.T)) xc_new = np.zeros(xc_old.shape) lc_new = np.zeros(lc_old.shape) for k in range(m): for i in range(n): for j in range(n): - if wc[i, j] != 0: # avoid divisions by extremely small values in wc_approx + if wc[i, j] != 0: # avoid divisions by zero in wc_approx xc_new[i, k] += wc[i, j] * lc_old[k, k] * xc_old[j, k] / wc_approx[i, j] lc_new[k, k] += wc[i, j] * xc_old[i, k] * xc_old[j, k] / wc_approx[i, j] xc_new[i, k] *= (2 * alpha * xc_old[i, k]) xc_new[i, k] += (1 - alpha) * yc[i, k] lc_new[k, k] *= (alpha * lc_old[k, k]) lc_new[k, k] += (1 - alpha) * sum(yc[:, k]) xc_new = normalize_cols(xc_new) lc_new /= np.trace(lc_new) gamma[it] = alpha * kl_divergence(wc, np.linalg.multi_dot([xc_new, lc_new, xc_new.T])) + (1 - alpha) * kl_divergence(yc, xc_new.dot(lc_new)) - # print("iter: {}, gamma: {}".format(it, gamma[it])) - if it == 0: gamma_min = gamma[it] else: - # print("d_gamma: {}, eps: {}, res: {}".format(abs((gamma[it] - gamma_min)) / gamma_min, eps, abs((gamma[it] - gamma_min)) / gamma_min < eps)) - if abs((gamma[it] - gamma_min)) / gamma_min < eps: - plt.plot(gamma[0:it]) - plt.ylabel(r"$\gamma$") - plt.xlabel("# iter") - plt.show() + if is_plot: + plt.plot(gamma[0:it]) + plt.ylabel(r"$\gamma$") + plt.xlabel("# iter") + plt.show() print("====> no iter: {}".format(it)) return xc_res, lc_res if gamma[it] < gamma[it - 1]: gamma_min = gamma[it] xc_res = xc_new lc_res = lc_new + else: + if is_plot: + plt.plot(gamma[0:it]) + plt.ylabel(r"$\gamma$") + plt.xlabel("# iter") + plt.show() + print("====> no iter: {}".format(it)) + return xc_res, lc_res + xc_old = xc_new lc_old = lc_new - raise Exception('Maximum iteration number reached: {}'.format(max_iter)) def read_edge_list(filename, weighted=False): - # TODO: I think the adjacency matrices come in unordered -> almost always 50 percent cluster association. print("Reading edge file {}".format(filename)) idmap = set() edge_cache = {} with open(filename) as f: for line in f: if weighted: u, v, w = [int(x) for x in line.strip().split()] else: tmp = [int(x) for x in line.strip().split()] u, v, w = tmp[0], tmp[1], 1.0 edge_cache[(u, v)] = w idmap.add(u) idmap.add(v) idmap = list(idmap) idmap_inv = {nid: i for i, nid in enumerate(idmap)} n = len(idmap) adj_mat = np.zeros((n, n)) for (u, v), w in edge_cache.items(): adj_mat[idmap_inv[u], idmap_inv[v]] = w adj_mat += adj_mat.T return idmap, idmap_inv, adj_mat -def alg(net_path, alpha, t_steps, n, m): - # FacetNet with fixed number of communities and individuals - - xc = np.random.rand(n, m) - xc = normalize_cols(xc) - lc = np.diag(np.random.rand(m)) - lc /= np.trace(lc) - - for t in range(t_steps): - idmap, idmap_inv, adj_mat = read_edge_list(net_path + "/%d.edgelist" % t, weighted=False) - - # TODO: this similarity calculation is experiment specific and must be done outside this function - # wc = similarity(adj_mat, weighted=True) - wc = normalize_rows(adj_mat) - - xc_old = xc - lc_old = lc - xc, lc = param_update(xc, lc, wc, alpha) - yc = xc.dot(lc) - - dc_inv = np.zeros(n) - for i in range(n): - dc_inv[i] = 1 / np.sum(yc[i, :]) - dc_inv = np.matrix(np.diag(dc_inv)) - - print(dc_inv) - print(yc) - - soft_comm = dc_inv.dot(yc) - - print("time:", t) - - print("soft_comm") - print(soft_comm) - - print("community net") - print(np.linalg.multi_dot([lc, xc.T, soft_comm])) - - print("evolution net") - print(np.linalg.multi_dot([lc_old, xc_old.T, soft_comm])) - - print("Activity (dc)") - print(np.diag(dc_inv)) +def facetnet_evolution(net_path, alpha, t_steps, m): + """ + Applies facetNet algorithm on an entire set of consecutively named edgelist.X files. Can deal with appearing and + disappearing number of nodes. Establishes soft community memberships to an a priori given number of communities (m).. - df = pd.DataFrame(data=soft_comm, columns=idmap) - df.to_csv("{}/soft_comm_nw{}".format(net_path, t)) - - -def alg_extended(net_path, alpha, t_steps, m): - # FacetNet with for variable number of nodes, a solution for different community numbers is not implemented + :param net_path: Path to edge-lists + :param alpha: Alpha cost weight + :param t_steps: Total number of time-steps to be computed (must be at least # of edge-list files) + :param m: Number of (a priori) communities + :return: No return value, saves soft_comm as soft_comm_alpha0.X_nwX.csv files in net_path folder + """ clusters = [] for i in range(m): clusters.append("cluster_{}".format(i)) for t in range(t_steps): print("time-step:", t) idmap, idmap_inv, adj_mat = read_edge_list(net_path + "/%d.edgelist" % t, weighted=True) n = len(idmap) - # wc = similarity(adj_mat, weighted=False) wc = normalize_rows(adj_mat) print("Calculate network representation ...") if t == 0: xc = np.random.rand(n, m) - xc = normalize_cols(xc) + xc = normalize_cols(xc) # TODO: Zero division happens here (step 7 row 43 in wc) lc = np.diag(np.random.rand(m)) lc /= np.trace(lc) yc = xc.dot(lc) xc, lc, = param_update(yc, wc, 1.0) else: reserved_rows = [idmap_inv0[x] for x in idmap0 if x in idmap] num_new = len(set(idmap) - set(idmap0)) num_old = len(reserved_rows) xc = xc[reserved_rows, :] xc = normalize_cols(xc) xc *= num_old / (num_old + num_new) if num_new != 0: xc = np.pad(xc, ((0, num_new), (0, 0)), mode='constant', constant_values=(1 / num_new, 1 / num_new)) yc = xc.dot(lc) xc, lc, = param_update(yc, wc, alpha) yc = xc.dot(lc) dc_inv = np.zeros(n) for i in range(n): dc_inv[i] = 1 / np.sum(yc[i, :]) dc_inv = np.matrix(np.diag(dc_inv)) soft_comm = dc_inv.dot(yc) df = pd.DataFrame(data=soft_comm, columns=clusters) df.insert(0, column="id", value=idmap) df.to_csv("{}/soft_comm_alpha{}_nw{}.csv".format(net_path, alpha, t), index=False) idmap0 = idmap idmap_inv0 = idmap_inv -def exp1(): - # Experiment with network stated in 4.1.2 - t_steps = 15 - from synthetic import generate_evolution_exp1 - print("generating synthetic graph") - generate_evolution_exp1("test_data/", tsteps=t_steps) - print("start the algorithm") - alpha = 0.5 - n, m = 128, 4 - np.random.seed(0) - alg("test_data/", alpha, t_steps, n, m) - - -def exp2(): - # Experiment with adding and removing nodes - t_steps = 15 - from synthetic import generate_evolution_exp2 - print("generating synthetic graph") - generate_evolution_exp2("./data/syntetic2/", tsteps=t_steps) - print("start the algorithm") - alpha = 0.5 - np.random.seed(0) - alg_extended("./data/syntetic2/", alpha, t_steps, 4) - - -def exp3(): - # Experiment with network stated in 4.1.2, adding weight - t_steps = 15 - from synthetic import generate_evolution_exp3 - print("generating synthetic graph") - generate_evolution_exp3("./data/syntetic3/", tsteps=t_steps) - print("start the algorithm") - alpha = 0.9 - n = 128 - m = 4 - np.random.seed(0) - alg("./data/syntetic3/", alpha, t_steps, n, m) - - if __name__ == "__main__": parser = argparse.ArgumentParser() parser.add_argument("t_steps", help="number of networks", type=int) - parser.add_argument("nw_folder", help="folder holding the network data", type=str) - parser.add_argument("alpha", help="alpha constant from gamma = CS * alpha + CT * (1 - alpha)", type=float) - parser.add_argument("m_cluster", help="number (m) of clusters", type=int) + parser.add_argument("alpha", help="alpha cost weight (0.0, 1.0] for gamma = CS*alpha + CT*(1-alpha)", type=float) + parser.add_argument("m_cluster", help="number of clusters", type=int) + parser.add_argument("nw_folder", help="folder holding the network data in edgelist.X form", type=str) args = parser.parse_args() - if args.alpha < 0.0 or args.alpha > 1.0: - print("Alpha value must be in [0.0, 1.0], terminate script") + if args.alpha <= 0.0 or args.alpha > 1.0: + print("Alpha value must be in (0.0, 1.0], terminate script") sys.exit(1) print("Analyzing networks from: {}".format(args.nw_folder)) print("# time-steps: {}, alpha: {}, # clusters: {}".format(args.t_steps, args.alpha, args.m_cluster)) - alg_extended(args.nw_folder, args.alpha, args.t_steps, args.m_cluster) - - print("soft_communities written to: {}".format(args.nw_folder)) - - # path = "/home/matthias/Documents/CodeDataRaph/data_Danielle/24h_networks/col4" - # exp_raph() - - # print("Experiment with network stated in 4.1.2") - # exp1() - - # print("Experiment with adding and removing nodes") - # exp2() + facetnet_evolution(args.nw_folder, args.alpha, args.t_steps, args.m_cluster) - # print("Experiment with network stated in 4.1.2, adding weight") - # exp3() + print("soft_communities written to: {}".format(args.nw_folder)) \ No newline at end of file