minor changes

Author: joshloyal

examples/multinet.py

@@ -14,9 +14,9 @@
 
 
 Y, X, lmbda, delta, probas, dists, z = correlated_dynamic_multilayer_network(
-    n_nodes=100, n_layers=4, n_time_steps=10, n_features=2,
-    center=0.75, tau=[0.25, 0.5], include_delta=True,
-    rho=0., rho_t=0.8, sigma=0.01,
+    n_nodes=100, n_layers=5, n_time_steps=10, n_features=2,
+    center=1, tau=[0.5, 0.75], include_delta=True,
+    rho=0., rho_t=0.4, sigma=0.01,
     random_state=1)
 print(lmbda)
 
@@ -35,7 +35,7 @@
 for t in range(n_time_steps):
     print('t = {}'.format(t))
     model = MultilayerNetworkLSM(
-                max_iter=50, n_features=2,
+                max_iter=500, n_features=2,
                 init_covariance_type='full',
                 tol=1e-2, n_init=10,
                 stopping_criteria='loglik',
@@ -55,7 +55,7 @@
 print(calculate_auc(Y, probas_pred, test_indices=test_indices))
 
 model_joint = DynamicMultilayerNetworkLSM(
-        max_iter=50, n_features=2,
+        max_iter=500, n_features=2,
         init_covariance_type='full',
         approx_type='structured',
         tol=1e-2, n_init=10,

examples/svi_mf.py

@@ -17,10 +17,10 @@
 sigma = 0.05
 
 Y, X, lmbda, delta, probas, dists, z = correlated_dynamic_multilayer_network(
-    n_nodes=100, n_layers=4, n_time_steps=10, n_features=2,
-    tau=np.array([0.5, 0.75]),
-    rho=0., rho_t=0.8, sigma=sigma,
-    random_state=3)
+    n_nodes=100, n_layers=5, n_time_steps=20,
+    n_features=2, tau=np.array([0.5, 0.75]),
+    center=1, rho_t=0.4, sigma=0.05,
+    random_state=83529)
 
 n_layers, n_time_steps, n_nodes, _ = Y.shape
 n_dyads = 0.5 * n_layers * n_time_steps * n_nodes * (n_nodes - 1)
@@ -29,17 +29,17 @@
 
 Y_train, test_indices = train_test_split(Y, test_size=0.2, random_state=23)
 
-#model_mf = DynamicMultilayerNetworkLSM(
-#                max_iter=50, n_features=2,
-#                init_covariance_type='spherical',
-#                approx_type='mean_field',
-#                tol=1e-2, n_init=10,
-#                stopping_criteria=None,
-#                n_jobs=-1,
-#                random_state=123)
-#
-#callback = TestMetricsCallback(Y=Y, probas=probas, test_indices=test_indices)
-#model_mf.fit(Y_train, callback=callback)
+model_mf = DynamicMultilayerNetworkLSM(
+                max_iter=50, n_features=2,
+                init_covariance_type='full',
+                approx_type='mean_field',
+                tol=1e-2, n_init=10,
+                stopping_criteria=None,
+                n_jobs=-1,
+                random_state=123)
+
+callback = TestMetricsCallback(Y=Y, probas=probas, test_indices=test_indices)
+model_mf.fit(Y_train, callback=callback)
 #print('MF Bias(sigma_sq): ', np.abs(model_mf.sigma_sq_ - sigma**2))
 #
 ## forecast accuracy
@@ -55,7 +55,11 @@
 #print('forecast AUC:', roc_auc_score(y_true, y_pred))
 #print('forecast correlation:', pearsonr(y_true_probas, y_pred)[0])
 #
-#plt.plot(model_mf.callback_.times_, model_mf.callback_.correlations_, label='mean_field')
+plt.plot(model_mf.callback_.times_, model_mf.callback_.correlations_,
+         linewidth=2, linestyle='dashed', label='mean_field')
+for i in range(10):
+    plt.plot(model_mf._callbacks[i].times_, model_mf._callbacks[i].correlations_,
+             alpha=0.25, c='black', linestyle='dashed')
 
 model_svi = DynamicMultilayerNetworkLSM(
                 max_iter=50, n_features=2,
@@ -70,13 +74,24 @@
 model_svi.fit(Y_train, callback=callback)
 #print('SVI Bias(sigma_sq): ', np.abs(model_svi.sigma_sq_ - sigma**2))
 #
-#plt.plot(model_svi.callback_.times_, model_svi.callback_.correlations_, label='structured')
+plt.plot(model_svi.callback_.times_, model_svi.callback_.correlations_, linewidth=2, label='structured')
+for i in range(10):
+    plt.plot(model_svi._callbacks[i].times_, model_svi._callbacks[i].correlations_,
+             alpha=0.25, c='black')
+
 #
 #plt.legend()
-#plt.xlabel('Wallclock Time [s]')
-#plt.ylabel('Pearson Correlation')
-#plt.show()
-#
+plt.xlabel('Wallclock Time [s]')
+plt.ylabel('Pearson Correlation')
+plt.show()
+
+
+plt.plot(model_mf.callback_.correlations_, label='Structured VI')
+plt.plot(model_svi.callback_.correlations_, label='Mean Field')
+plt.xlabel('# of Iterations')
+plt.ylabel('Pearson Correlation')
+plt.show()
+
 ## forecast accuracy
 #y_true = []
 #y_true_probas = []

multidynet/datasets/synthetic.py

@@ -176,7 +176,7 @@ def correlated_dynamic_multilayer_network(n_nodes=100, n_layers=4, n_time_steps=
             X[t] = X[t-1] + errors[..., t-1]
             X[t] -= np.mean(X[t], axis=0)
 
-        # sample assortativity parameters from a U(-2, 2)
+        # sample assortativity parameters from a U(-1, 2) (average is 0.5)
         lmbda = np.zeros((n_layers, n_features))
         lmbda[0] = rng.choice([-1, 1], size=n_features)
         lmbda[1:] = rng.uniform(

multidynet/metrics.py

@@ -7,6 +7,52 @@
 from sklearn.metrics import roc_auc_score
 
 
+def calculate_metric_single(Y_true, Y_pred, metric, test_indices=None):
+    n_time_steps, n_nodes, _ = Y_true.shape
+    indices = np.tril_indices_from(Y_true[0], k=-1)
+
+    y_true = []
+    y_pred = []
+    for t in range(n_time_steps):
+        y_true_vec = Y_true[t][indices]
+        y_pred_vec = Y_pred[t][indices]
+
+        if test_indices is None:
+            subset = y_true_vec != -1.0
+        else:
+            subset = test_indices[t]
+        y_true.extend(y_true_vec[subset])
+        y_pred.extend(y_pred_vec[subset])
+
+    return metric(y_true,  y_pred)
+
+
+def calculate_metric(Y_true, Y_pred, metric, test_indices=None):
+    if Y_true.ndim == 3:
+        return calculate_metric_single(
+                Y_true, Y_pred, metric, test_indices=test_indices)
+
+    n_layers, n_time_steps, n_nodes, _ = Y_true.shape
+    indices = np.tril_indices_from(Y_true[0, 0], k=-1)
+
+    y_true = []
+    y_pred = []
+    for k in range(n_layers):
+        for t in range(n_time_steps):
+            y_true_vec = Y_true[k, t][indices]
+            y_pred_vec = Y_pred[k, t][indices]
+
+            if test_indices is None:
+                subset = y_true_vec != -1.0
+            else:
+                subset = test_indices[k][t]
+
+            y_true.extend(y_true_vec[subset])
+            y_pred.extend(y_pred_vec[subset])
+
+    return metric(y_true, y_pred)
+
+
 def calculate_auc_single(Y_true, Y_pred, test_indices=None):
     n_time_steps, n_nodes, _ = Y_true.shape
     indices = np.tril_indices_from(Y_true[0], k=-1)

multidynet/multidynet.py

@@ -7,8 +7,9 @@
 from scipy.special import logit, gammainc, expit, logsumexp
 from sklearn.exceptions import ConvergenceWarning
 from sklearn.utils import check_array, check_random_state
+from scipy.optimize import linear_sum_assignment
 from sklearn.linear_model import LogisticRegression
-from sklearn.metrics import log_loss
+from sklearn.metrics import log_loss, f1_score
 from tqdm import tqdm
 
 from .omega import update_omega
@@ -18,7 +19,7 @@
 from .variances import update_tau_sq, update_sigma_sq, update_X0_precision, update_diag_tau_sq
 from .variances import update_tau_sq_delta, update_sigma_sq_delta
 from .log_likelihood import log_likelihood, pointwise_log_likelihood
-from .metrics import calculate_auc
+from .metrics import calculate_auc, calculate_metric
 from .sample_lds import sample_gssm
 from .model_selection import dynamic_multilayer_adjacency_to_vec
 
@@ -60,6 +61,22 @@ def __init__(self, omega, X, X_sigma, X_cross_cov,
         self.callback_ = callback
 
 
+def find_permutation(U, U_ref):
+    _, n_features = U.shape
+    C = U_ref.T @ U
+    perm = linear_sum_assignment(np.maximum(C, -C), maximize=True)[1]
+    sign = np.sign(C[np.arange(n_features), perm])
+    return sign * U[:, perm]
+
+
+def smooth_positions(U):
+    n_time_steps, _, _ = U.shape
+    for t in range(1, n_time_steps):
+        U[t] = find_permutation(U[t], U[t-1])
+
+    return U
+
+
 def sample_socialities(model, size=500, random_state=None):
     samples = model.sample(size=size, random_state=random_state)
     deltas = samples['delta']
@@ -556,6 +573,9 @@ def fit(self, Y, callback=None, n_samples=None):
 
         self._set_parameters(best_model)
 
+        # save all callbacks
+        self._callbacks = [models[i].callback_ for i in range(len(models))]
+
         # calculate dyad-probabilities
         self.probas_ = calculate_probabilities(
             self.X_, self.lambda_, self.delta_)
@@ -659,6 +679,9 @@ def _set_parameters(self, model):
         self.X_sigma_ = model.X_sigma_
         self.X_cross_cov_ = model.X_cross_cov_
 
+        # match signed-permutations across time
+        self.X_ = smooth_positions(self.X_)
+
         # transform to identifiable parameterization
         if self.X_ is not None:
             self.Z_ = (self.X_ -