From 3e819e8816f2099dc17871d534de06d432e8a744 Mon Sep 17 00:00:00 2001 From: nicoval Date: Fri, 5 Mar 2021 15:25:18 +0100 Subject: [PATCH 1/4] v2.0.2: solve standard input initial_guess problem in solve_tool --- src/NEMtropy/graph_classes.py | 8 ++++---- tests/test_solve_tool.py | 2 -- 2 files changed, 4 insertions(+), 6 deletions(-) diff --git a/src/NEMtropy/graph_classes.py b/src/NEMtropy/graph_classes.py index 77a642ea..6ddaf269 100644 --- a/src/NEMtropy/graph_classes.py +++ b/src/NEMtropy/graph_classes.py @@ -924,8 +924,8 @@ def _set_solved_problem_ecm(self, solution): def solve_tool( self, model, - method, - initial_guess=None, + method='newton', + initial_guess='random', adjacency="cm_exp", method_adjacency="newton", initial_guess_adjacency="random", @@ -2380,8 +2380,8 @@ def _set_solved_problem_crema_directed(self, solution): def solve_tool( self, model, - method, - initial_guess=None, + method='newton', + initial_guess='random', adjacency=None, method_adjacency='newton', initial_guess_adjacency="random", diff --git a/tests/test_solve_tool.py b/tests/test_solve_tool.py index bb431f82..0fe19df2 100644 --- a/tests/test_solve_tool.py +++ b/tests/test_solve_tool.py @@ -27,8 +27,6 @@ def test_dcm(self): g.solve_tool( model="dcm", - method="quasinewton", - initial_guess="uniform", max_steps=200, verbose=False, ) From c8c09b9b1f357f9ea00124df57b5ccafb610b430 Mon Sep 17 00:00:00 2001 From: nicoval Date: Fri, 5 Mar 2021 16:33:03 +0100 Subject: [PATCH 2/4] print solution error always in solve_tool --- src/NEMtropy/graph_classes.py | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/src/NEMtropy/graph_classes.py b/src/NEMtropy/graph_classes.py index 6ddaf269..01b374d1 100644 --- a/src/NEMtropy/graph_classes.py +++ b/src/NEMtropy/graph_classes.py @@ -778,7 +778,7 @@ def _initialize_problem(self, model, method): if self.regularise == "eigenvalues": self.hessian_regulariser = sof.matrix_regulariser_function_eigen_based elif self.regularise == "identity": - self.hessian_regulariser = hessian_regulariser_function + self.hessian_regulariser = sof.hessian_regulariser_function def _solve_problem_crema_undirected( self, @@ -1020,8 +1020,7 @@ def solve_tool( eps=eps, ) self._solution_error() - if verbose: - print("\nmin eig = {}".format(self.error)) + print("\nsolution error = {}".format(self.error)) def ensemble_sampler(self, n, cpu_n=1, output_dir="sample/", seed=42): """The function sample a given number of graphs in the ensemble @@ -2476,8 +2475,7 @@ def solve_tool( eps=eps, ) self._solution_error() - if verbose: - print("\nmin eig = {}".format(self.error)) + print("\nsolution error = {}".format(self.error)) def ensemble_sampler(self, n, cpu_n=1, output_dir="sample/", seed=42): """The function sample a given number of graphs in the ensemble From 0761568a9b38591a43717e6a54ed1473496940cc Mon Sep 17 00:00:00 2001 From: nicoval Date: Fri, 26 Mar 2021 19:41:40 +0100 Subject: [PATCH 3/4] fixed default seed setting --- src/NEMtropy/graph_classes.py | 4 +- tests/test_ensemble_seed.py | 96 +++++++++++++++++++++++++++++++++++ 2 files changed, 98 insertions(+), 2 deletions(-) create mode 100644 tests/test_ensemble_seed.py diff --git a/src/NEMtropy/graph_classes.py b/src/NEMtropy/graph_classes.py index 01b374d1..e8a4d3a8 100644 --- a/src/NEMtropy/graph_classes.py +++ b/src/NEMtropy/graph_classes.py @@ -1022,7 +1022,7 @@ def solve_tool( self._solution_error() print("\nsolution error = {}".format(self.error)) - def ensemble_sampler(self, n, cpu_n=1, output_dir="sample/", seed=42): + def ensemble_sampler(self, n, cpu_n=1, output_dir="sample/", seed=None): """The function sample a given number of graphs in the ensemble generated from the last model solved. Each grpah is an edgelist written in the output directory as `.txt` file. @@ -2477,7 +2477,7 @@ def solve_tool( self._solution_error() print("\nsolution error = {}".format(self.error)) - def ensemble_sampler(self, n, cpu_n=1, output_dir="sample/", seed=42): + def ensemble_sampler(self, n, cpu_n=1, output_dir="sample/", seed=None): """The function sample a given number of graphs in the ensemble generated from the last model solved. Each grpah is an edgelist written in the output directory as `.txt` file. diff --git a/tests/test_ensemble_seed.py b/tests/test_ensemble_seed.py new file mode 100644 index 00000000..ce2b7247 --- /dev/null +++ b/tests/test_ensemble_seed.py @@ -0,0 +1,96 @@ +"""This test checks the default setting for the seed is working. +When the seed is not manually set, it should be ranodmly chosen. +""" + +import sys +import os +sys.path.append("../") +import NEMtropy.graph_classes as sample +import NEMtropy.matrix_generator as mg +import numpy as np +import unittest # test tool +import random +import networkx as nx + + +class MyTest(unittest.TestCase): + + + def setUp(self): + pass + + + def test_0(self): + """test with 3 classes of cardinality 1 + and no zero degrees + """ + """ + A = np.array( + [ + [0, 1, 1, 0], + [1, 0, 0, 1], + [1, 0, 0, 0], + [0, 1, 0, 0], + ] + ) + e = [(0,1), (0,2), (1,3)] + d = [1,1,2,2] + print(e) + print(d) + """ + N, seed = (10, 42) + A = mg.random_binary_matrix_generator_dense(N, sym=True, seed=seed) + # number of copies to generate + + g = sample.UndirectedGraph(A) + + g._solve_problem( + model="cm", + method="fixed-point", + max_steps=100, + verbose=False, + linsearch=True, + initial_guess="uniform", + ) + + x = g.x + # g._solution_error() + err = g.error + + # print('\ntest 5: error = {}'.format(g.error)) + n = 10 + output_dir = "sample/" + # random.seed(100) + g_list = [] + for i in range(n): + g.ensemble_sampler(n=1, output_dir=output_dir) + g_list.append(np.loadtxt("sample/0.txt")) + + appo = True + old = g_list[0] + for i in range(1, n): + appo = appo*np.all(old == g_list[i]) + + # debug + """ + print('original dseq',d) + print('original dseq sum ',g.dseq.sum()) + print('ensemble average dseq', d_emp) + print('ensemble dseq sum ',np.array([d_emp[key] for key in d_emp.keys()]).sum()) + print(d_diff) + print('empirical error', ensemble_error) + print('theoretical error', err) + """ + + + l = os.listdir(output_dir) + for f in l: + os.remove(output_dir + f) + os.rmdir(output_dir) + + # test result + self.assertTrue(not appo) + + +if __name__ == "__main__": + unittest.main() From a1bb7a632f14da8b7c6818495031aec0852c5f61 Mon Sep 17 00:00:00 2001 From: nicoval Date: Tue, 30 Mar 2021 17:20:31 +0200 Subject: [PATCH 4/4] added model_loglikelihood method and test --- src/NEMtropy/graph_classes.py | 93 +++++++++++------- tests/test_model_loglikelihood.py | 153 ++++++++++++++++++++++++++++++ 2 files changed, 211 insertions(+), 35 deletions(-) create mode 100644 tests/test_model_loglikelihood.py diff --git a/src/NEMtropy/graph_classes.py b/src/NEMtropy/graph_classes.py index e8a4d3a8..b85b7add 100644 --- a/src/NEMtropy/graph_classes.py +++ b/src/NEMtropy/graph_classes.py @@ -88,6 +88,7 @@ def __init__( self.relative_error_strength = None self.full_return = False self.last_model = None + self.solution_array # function self.args = None @@ -334,16 +335,16 @@ def _solve_problem( def _set_solved_problem_cm(self, solution): if self.full_return: - self.r_xy = solution[0] + self.solution_array = solution[0] self.comput_time = solution[1] self.n_steps = solution[2] self.norm_seq = solution[3] self.diff_seq = solution[4] self.alfa_seq = solution[5] else: - self.r_xy = solution + self.solution_array = solution - self.r_x = self.r_xy + self.r_x = self.solution_array if self.last_model == "cm": self.x = self.r_x[self.r_invert_dseq] elif self.last_model == "cm_exp": @@ -903,23 +904,25 @@ def _set_solved_problem_crema_undirected(self, solution): else: self.beta = solution + self.solution_array = self.beta + def _set_solved_problem_ecm(self, solution): if self.full_return: - self.r_xy = solution[0] + self.solution_array = solution[0] self.comput_time = solution[1] self.n_steps = solution[2] self.norm_seq = solution[3] self.diff_seq = solution[4] self.alfa_seq = solution[5] else: - self.r_xy = solution + self.solution_array = solution if self.last_model == "ecm": - self.x = self.r_xy[: self.n_nodes] - self.y = self.r_xy[self.n_nodes:] + self.x = self.solution_array[: self.n_nodes] + self.y = self.solution_array[self.n_nodes:] elif self.last_model == "ecm_exp": - self.x = np.exp(-self.r_xy[:self.n_nodes]) - self.y = np.exp(-self.r_xy[self.n_nodes:]) + self.x = np.exp(-self.solution_array[:self.n_nodes]) + self.y = np.exp(-self.solution_array[self.n_nodes:]) def solve_tool( self, @@ -1120,6 +1123,12 @@ def ensemble_sampler(self, n, cpu_n=1, output_dir="sample/", seed=None): else: raise ValueError("insert a model") + def model_loglikelihood(self): + """Returns the loglikelihood of the solution of last model executed. + """ + return self.step_fun(self.solution_array) + + class DirectedGraph: """Directed graph instance can be initialised with @@ -1188,7 +1197,7 @@ def __init__( # reduced solutions self.r_x = None self.r_y = None - self.r_xy = None + self.solution_array = None # Problem (reduced) residuals self.residuals = None self.final_result = None @@ -1478,17 +1487,17 @@ def _solve_problem( def _set_solved_problem_dcm(self, solution): if self.full_return: - self.r_xy = solution[0] + self.solution_array = solution[0] self.comput_time = solution[1] self.n_steps = solution[2] self.norm_seq = solution[3] self.diff_seq = solution[4] self.alfa_seq = solution[5] else: - self.r_xy = solution + self.solution_array = solution - self.r_x = self.r_xy[: self.rnz_n_out] - self.r_y = self.r_xy[self.rnz_n_out:] + self.r_x = self.solution_array[: self.rnz_n_out] + self.r_y = self.solution_array[self.rnz_n_out:] self.x = self.r_x[self.r_invert_dseq] self.y = self.r_y[self.r_invert_dseq] @@ -1496,7 +1505,7 @@ def _set_solved_problem_dcm(self, solution): def _set_solved_problem_dcm_exp(self, solution): if self.full_return: # conversion from theta to x - self.r_xy = np.exp(-solution[0]) + self.solution_array = np.exp(-solution[0]) self.comput_time = solution[1] self.n_steps = solution[2] self.norm_seq = solution[3] @@ -1504,34 +1513,34 @@ def _set_solved_problem_dcm_exp(self, solution): self.alfa_seq = solution[5] else: # conversion from theta to x - self.r_xy = np.exp(-solution) + self.solution_array = np.exp(-solution) - self.r_x = self.r_xy[: self.rnz_n_out] - self.r_y = self.r_xy[self.rnz_n_out:] + self.r_x = self.solution_array[: self.rnz_n_out] + self.r_y = self.solution_array[self.rnz_n_out:] self.x = self.r_x[self.r_invert_dseq] self.y = self.r_y[self.r_invert_dseq] def _set_solved_problem_decm(self, solution): if self.full_return: - self.r_xy = solution[0] + self.solution_array = solution[0] self.comput_time = solution[1] self.n_steps = solution[2] self.norm_seq = solution[3] self.diff_seq = solution[4] self.alfa_seq = solution[5] else: - self.r_xy = solution + self.solution_array = solution - self.x = self.r_xy[: self.n_nodes] - self.y = self.r_xy[self.n_nodes: 2 * self.n_nodes] - self.b_out = self.r_xy[2 * self.n_nodes: 3 * self.n_nodes] - self.b_in = self.r_xy[3 * self.n_nodes:] + self.x = self.solution_array[: self.n_nodes] + self.y = self.solution_array[self.n_nodes: 2 * self.n_nodes] + self.b_out = self.solution_array[2 * self.n_nodes: 3 * self.n_nodes] + self.b_in = self.solution_array[3 * self.n_nodes:] def _set_solved_problem_decm_exp(self, solution): if self.full_return: # conversion from theta to x - self.r_xy = np.exp(-solution[0]) + self.solution_array = np.exp(-solution[0]) self.comput_time = solution[1] self.n_steps = solution[2] self.norm_seq = solution[3] @@ -1539,12 +1548,12 @@ def _set_solved_problem_decm_exp(self, solution): self.alfa_seq = solution[5] else: # conversion from theta to x - self.r_xy = np.exp(-solution) + self.solution_array = np.exp(-solution) - self.x = self.r_xy[: self.n_nodes] - self.y = self.r_xy[self.n_nodes: 2 * self.n_nodes] - self.b_out = self.r_xy[2 * self.n_nodes: 3 * self.n_nodes] - self.b_in = self.r_xy[3 * self.n_nodes:] + self.x = self.solution_array[: self.n_nodes] + self.y = self.solution_array[self.n_nodes: 2 * self.n_nodes] + self.b_out = self.solution_array[2 * self.n_nodes: 3 * self.n_nodes] + self.b_in = self.solution_array[3 * self.n_nodes:] def _set_solved_problem(self, solution): model = self.last_model @@ -2376,6 +2385,8 @@ def _set_solved_problem_crema_directed(self, solution): self.b_out = solution[: self.n_nodes] self.b_in = solution[self.n_nodes:] + self.solution_array = solution + def solve_tool( self, model, @@ -2579,6 +2590,11 @@ def ensemble_sampler(self, n, cpu_n=1, output_dir="sample/", seed=None): else: raise ValueError("insert a model") + def model_loglikelihood(self): + """Returns the loglikelihood of the solution of last model executed. + """ + return self.step_fun(self.solution_array) + class BipartiteGraph: """Bipartite Graph class for undirected binary bipartite networks. @@ -2639,7 +2655,7 @@ def __init__(self, biadjacency=None, adjacency_list=None, edgelist=None, degree_ self.y = None self.r_x = None self.r_y = None - self.r_xy = None + self.solution_array = None self.dict_x = None self.dict_y = None self.theta_x = None @@ -2681,6 +2697,7 @@ def __init__(self, biadjacency=None, adjacency_list=None, edgelist=None, degree_ self.full_rows_num = None self.full_rows_num = None self.solution_converged = None + self.solution_array = None def _initialize_graph(self, biadjacency=None, adjacency_list=None, edgelist=None, degree_sequences=None): """ @@ -3137,15 +3154,15 @@ def _set_solved_problem(self, solution): self.r_theta_xy = solution self.r_theta_x = self.r_theta_xy[:self.r_n_rows] self.r_theta_y = self.r_theta_xy[self.r_n_rows:] - self.r_xy = np.exp(- self.r_theta_xy) + self.solution_array = np.exp(- self.r_theta_xy) self.r_x = np.exp(- self.r_theta_x) self.r_y = np.exp(- self.r_theta_y) self.theta_x[self.nonfixed_rows] = self.r_theta_x[self.r_invert_rows_deg] self.theta_y[self.nonfixed_cols] = self.r_theta_y[self.r_invert_cols_deg] else: - self.r_xy = solution - self.r_x = self.r_xy[:self.r_n_rows] - self.r_y = self.r_xy[self.r_n_rows:] + self.solution_array = solution + self.r_x = self.solution_array[:self.r_n_rows] + self.r_y = self.solution_array[self.r_n_rows:] if self.x is None: self.x = np.zeros(self.n_rows) if self.y is None: @@ -3686,3 +3703,9 @@ def clean_edges(self): self.rows_deg = None self.cols_deg = None self.is_initialized = False + + + def model_loglikelihood(self): + """Returns the loglikelihood of the solution of last model executed. + """ + return self.step_fun(self.solution_array) diff --git a/tests/test_model_loglikelihood.py b/tests/test_model_loglikelihood.py new file mode 100644 index 00000000..8287df77 --- /dev/null +++ b/tests/test_model_loglikelihood.py @@ -0,0 +1,153 @@ +import sys + +sys.path.append("../") +import NEMtropy.graph_classes as sample +import numpy as np +import unittest # test tool +import NEMtropy.matrix_generator as mg +import NEMtropy.models_functions as mof + + +class MyTest(unittest.TestCase): + def setUp(self): + pass + + def test_dcm(self): + """test with 3 classes of cardinality 1 + and no zero degrees + """ + A = np.array( + [ + [0, 1, 1], + [1, 0, 0], + [0, 1, 0], + ] + ) + + g = sample.DirectedGraph(A) + + g.solve_tool( + model="dcm", + max_steps=200, + verbose=False, + ) + l = -mof.loglikelihood_dcm(g.solution_array, g.args) + + # g._solution_error() + # debug + # print(g.r_dseq_out) + # print(g.r_dseq_in) + # print(g.rnz_dseq_out) + # print(g.rnz_dseq_in) + # print('\ntest 0: error = {}'.format(g.error)) + + # test result + self.assertTrue(l == g.model_loglikelihood()) + + def test_dcm_new(self): + """test with 3 classes of cardinality 1 + and no zero degrees + """ + A = np.array( + [ + [0, 1, 1], + [1, 0, 0], + [0, 1, 0], + ] + ) + + g = sample.DirectedGraph(A) + + g.solve_tool( + model="dcm_exp", + method="quasinewton", + initial_guess="uniform", + max_steps=200, + verbose=False, + ) + l = -mof.loglikelihood_dcm_exp(g.solution_array, g.args) + + # g._solution_error() + # debug + # print(g.r_dseq_out) + # print(g.r_dseq_in) + # print(g.rnz_dseq_out) + # print(g.rnz_dseq_in) + # print('\ntest 0: error = {}'.format(g.error)) + # print(l) + # print(g.model_loglikelihood()) + + # test result + self.assertTrue(l == g.model_loglikelihood()) + + @unittest.skip("") + def test_decm(self): + """test with 3 classes of cardinality 1 + and no zero degrees + """ + + n, s = (4, 25) + + A = mg.random_weighted_matrix_generator_dense( + n, sup_ext=10, sym=False, seed=s, intweights=True + ) + + g = sample.DirectedGraph(A) + + g.solve_tool( + model="decm", + method="quasinewton", + initial_guess="uniform", + max_steps=200, + verbose=False, + ) + + l = -mof.loglikelihood_decm(g.solution_array, g.args) + + # g._solution_error() + # debug + # print(g.r_dseq_out) + # print(g.r_dseq_in) + # print(g.rnz_dseq_out) + # print(g.rnz_dseq_in) + # print('\ntest 0: error = {}'.format(g.error)) + + # test result + self.assertTrue(l == g.model_loglikelihood()) + + def test_decm_new(self): + """test with 3 classes of cardinality 1 + and no zero degrees + """ + + n, s = (10, 25) + + A = mg.random_weighted_matrix_generator_dense( + n, sup_ext=10, sym=False, seed=s, intweights=True + ) + + g = sample.DirectedGraph(A) + + g.solve_tool( + model="decm_exp", + method="quasinewton", + initial_guess="uniform", + max_steps=200, + verbose=False, + ) + l = -mof.loglikelihood_decm_exp(g.solution_array, g.args) + + # g._solution_error() + # debug + # print(g.r_dseq_out) + # print(g.r_dseq_in) + # print(g.rnz_dseq_out) + # print(g.rnz_dseq_in) + # print('\ntest 0: error = {}'.format(g.error)) + + # test result + self.assertTrue(l == g.model_loglikelihood()) + + +if __name__ == "__main__": + unittest.main()