diff --git a/multidynet/lmbdas_svi.pyx b/multidynet/lmbdas_svi.pyx new file mode 100644 index 0000000..d2302ad --- /dev/null +++ b/multidynet/lmbdas_svi.pyx @@ -0,0 +1,380 @@ +# encoding: utf-8 +# cython: language_level=3 +# cython: cdivision=True +# cython: boundscheck=False +# cython: wraparound=False +# cython: nonecheck=False +# cython: initializedcheck=False + +from libc.math cimport sqrt, ceil +from scipy.special.cython_special import expit + +import numpy as np +cimport numpy as np + +from .omega cimport update_omega_svi + + +ctypedef np.int64_t INT + + +cdef struct Index: + size_t i + size_t j + + +cdef inline Index get_row_col(int idx): + cdef Index index + + index.i = (ceil(sqrt(2 * (idx + 1) + 0.25) - 0.5)) + index.j = (idx - 0.5 * index.i * (index.i - 1)) + + return index + + +class LambdaParams(object): + def __init__(self, + lmbda, + lmbda_sigma, + lmbda0_eta, + lmbda_eta1, + lmbda_eta2): + self.lmbda = lmbda + self.lmbda_sigma = lmbda_sigma + self.lmbda0_eta = lmbda0_eta + self.lmbda_eta1 = lmbda_eta1 + self.lmbda_eta2 = lmbda_eta2 + + + +def sample_omegas(double[:, :, ::1] X, + double[:, :, :, ::1] X_sigma, + double[:, ::1] lmbda, + double[:, :, ::1] lmbda_sigma, + double[:, ::1] delta, + double[:, ::1] delta_sigma, + INT[:, :, ::1] pos_subsamples, + INT[:, :, ::1] neg_subsamples): + cdef size_t k, t, l, i, j + cdef int idx + cdef size_t n_time_steps = X.shape[0] + cdef size_t n_layers = lmbda.shape[0] + cdef size_t n_pos_subsamples = pos_subsamples.shape[2] + cdef size_t n_neg_subsamples = neg_subsamples.shape[2] + cdef Index index + + cdef np.ndarray[double, ndim=3, mode='c'] omega_pos = np.zeros( + (n_layers, n_time_steps, n_pos_subsamples), np.float64) + cdef np.ndarray[double, ndim=3, mode='c'] omega_neg = np.zeros( + (n_layers, n_time_steps, n_neg_subsamples), np.float64) + + for k in range(n_layers): + for t in range(n_time_steps): + + # loop over positive dyads + for l in range(n_pos_subsamples): + idx = pos_subsamples[k, t, l] + if idx < 0: + break + + index = get_row_col(idx) + i = index.i + j = index.j + + omega_pos[k, t, l] = update_omega_svi( + X[t, i], X_sigma[t, i], X[t, j], X_sigma[t, j], + lmbda[k], lmbda_sigma[k], delta[k, i], delta[k, j], + delta_sigma[k, i], delta_sigma[k, j]) + + # loop over positive dyads + for l in range(n_neg_subsamples): + idx = neg_subsamples[k, t, l] + if idx < 0: + break + + index = get_row_col(idx) + i = index.i + j = index.j + + omega_neg[k, t, l] = update_omega_svi( + X[t, i], X_sigma[t, i], X[t, j], X_sigma[t, j], + lmbda[k], lmbda_sigma[k], delta[k, i], delta[k, j], + delta_sigma[k, i], delta_sigma[k, j]) + + return omega_pos, omega_neg + + +cpdef double calculate_natural_parameters_reference(const double[:, :, :, ::1] Y, + double[:, :, ::1] X, + double[:, :, :, ::1] X_sigma, + double[:, :] delta, + double[:] lmbda, + double[:, :, ::1] omega_pos, + double[:, :, ::1] omega_neg, + INT[:, :, ::1] pos_subsamples, + INT[:, :, ::1] neg_subsamples, + double[:, ::1] pos_ratios, + double[:, ::1] neg_ratios, + int p): + cdef size_t t, i, l, q + cdef int idx + cdef size_t n_time_steps = Y.shape[1] + cdef size_t n_nodes = X.shape[1] + cdef size_t n_features = X.shape[2] + cdef size_t n_pos_subsamples = pos_subsamples.shape[2] + cdef size_t n_neg_subsamples = neg_subsamples.shape[2] + cdef double ratio + cdef Index index + cdef double eta = 0. + + for t in range(n_time_steps): + # loop over positive dyads + for l in range(n_pos_subsamples): + idx = pos_subsamples[0, t, l] + if idx < 0: + break + + index = get_row_col(idx) + i = index.i + j = index.j + ratio = pos_ratios[0, t] + + eta += ratio * ( + (Y[0, t, i, j] - 0.5 - omega_pos[0, t, l] * ( + delta[0, i] + delta[0, j])) * + X[t, i, p] * X[t, j, p]) + + for q in range(n_features): + if q != p: + eta -= ratio * omega_pos[0, t, l] * lmbda[q] * ( + (X_sigma[t, i, q, p] + X[t, i, q] * X[t, i, p]) * + (X_sigma[t, j, q, p] + X[t, j, q] * X[t, j, p])) + + # loop over negative dyads + for l in range(n_neg_subsamples): + idx = neg_subsamples[0, t, l] + if idx < 0: + break + + index = get_row_col(idx) + i = index.i + j = index.j + ratio = neg_ratios[0, t] + + eta += ratio * ( + (Y[0, t, i, j] - 0.5 - omega_neg[0, t, l] * ( + delta[0, i] + delta[0, j])) * + X[t, i, p] * X[t, j, p]) + + for q in range(n_features): + if q != p: + eta -= ratio * omega_neg[0, t, l] * lmbda[q] * ( + (X_sigma[t, i, q, p] + X[t, i, q] * X[t, i, p]) * + (X_sigma[t, j, q, p] + X[t, j, q] * X[t, j, p])) + + return 2 * eta + + +def calculate_natural_parameters(const double[:, :, :, ::1] Y, + double[:, :, ::1] X, + double[:, :, :, ::1] X_sigma, + double[:, :] delta, + double lmbda_var_prior, + double[:, :, ::1] omega_pos, + double[:, :, ::1] omega_neg, + INT[:, :, ::1] pos_subsamples, + INT[:, :, ::1] neg_subsamples, + double[:, ::1] pos_ratios, + double[:, ::1] neg_ratios, + int k): + cdef size_t t, i, j, l, p, q + cdef int idx + cdef size_t n_time_steps = Y.shape[1] + cdef size_t n_nodes = X.shape[1] + cdef size_t n_features = X.shape[2] + cdef size_t n_pos_subsamples = pos_subsamples.shape[2] + cdef size_t n_neg_subsamples = neg_subsamples.shape[2] + cdef double ratio + cdef Index index + + cdef np.ndarray[double, ndim=1, mode='c'] eta1 = np.zeros(n_features) + cdef np.ndarray[double, ndim=2, mode='c'] eta2 = np.zeros( + (n_features, n_features)) + + for t in range(n_time_steps): + # loop over positive dyads + for l in range(n_pos_subsamples): + idx = pos_subsamples[k, t, l] + if idx < 0: + break + + index = get_row_col(idx) + i = index.i + j = index.j + ratio = pos_ratios[k, t] + + for p in range(n_features): + eta1[p] += ratio * ( + (Y[k, t, i, j] - 0.5 - + omega_pos[k, t, l] * ( + delta[k, i] + delta[k, j])) * + X[t, i, p] * X[t, j, p]) + + for q in range(p + 1): + eta2[p, q] += ratio * omega_pos[k, t, l] * ( + (X_sigma[t, i, p, q] + X[t, i, p] * X[t, i, q]) * + (X_sigma[t, j, p, q] + X[t, j, p] * X[t, j, q])) + eta2[q, p] = eta2[p, q] + + # loop over negative dyads + for l in range(n_neg_subsamples): + idx = neg_subsamples[k, t, l] + if idx < 0: + break + + index = get_row_col(idx) + i = index.i + j = index.j + ratio = neg_ratios[k, t] + + for p in range(n_features): + eta1[p] += ratio * ( + (Y[k, t, i, j] - 0.5 - + omega_neg[k, t, l] * ( + delta[k, i] + delta[k, j])) * + X[t, i, p] * X[t, j, p]) + + for q in range(p + 1): + eta2[p, q] += ratio * omega_neg[k, t, l] * ( + (X_sigma[t, i, p, q] + X[t, i, p] * X[t, i, q]) * + (X_sigma[t, j, p, q] + X[t, j, p] * X[t, j, q])) + eta2[q, p] = eta2[p, q] + + eta2[np.diag_indices_from(eta2)] += (1. / lmbda_var_prior) + + return eta1, eta2 + + +def subsample_dyads(np.ndarray[double, ndim=4, mode='c'] Y, + int max_edges, + int n_neg_multiplier, + object rng): + cdef size_t k, t + cdef size_t n_layers = Y.shape[0] + cdef size_t n_time_steps = Y.shape[1] + cdef size_t n_nodes = Y.shape[2] + cdef size_t n_pos_samples, n_neg_samples, n_neg_subsamples + + cdef np.ndarray[double, ndim=1, mode='c'] y_vec + cdef np.ndarray[np.int64_t, ndim=1, mode='c'] indices + cdef np.ndarray[np.int64_t, ndim=3, mode='c'] pos_subsamples = np.full( + (n_layers, n_time_steps, max_edges), -1, dtype=np.int64) + cdef np.ndarray[double, ndim=2, mode='c'] pos_ratios = np.empty( + (n_layers, n_time_steps), dtype=np.float64) + cdef np.ndarray[np.int64_t, ndim=3, mode='c'] neg_subsamples = np.full( + (n_layers, n_time_steps, n_neg_multiplier * max_edges), -1, dtype=np.int64) + cdef np.ndarray[double, ndim=2, mode='c'] neg_ratios = np.empty( + (n_layers, n_time_steps), dtype=np.float64) + + for k in range(n_layers): + for t in range(n_time_steps): + y_vec = Y[k, t][np.tril_indices_from(Y[k, t], k=-1)] + + # XXX: positive dyads could be cached + indices = np.where(y_vec == 1)[0] + n_pos_samples = indices.shape[0] + if n_pos_samples > 0: + pos_subsamples[k, t, :indices.shape[0]] = indices + pos_ratios[k, t] = n_pos_samples / ( indices.shape[0]) + + # XXX: negative dyads could be cached + indices = np.where(y_vec == 0)[0] + n_neg_samples = indices.shape[0] + if n_neg_samples > 0: + n_neg_subsamples = min(n_neg_samples, n_neg_multiplier * n_pos_samples) + if n_neg_subsamples < n_neg_samples: + indices = rng.choice(indices, size=n_neg_subsamples) + neg_subsamples[k, t, :indices.shape[0]] = indices + neg_ratios[k, t] = n_neg_samples / ( indices.shape[0]) + + + return pos_subsamples, neg_subsamples, pos_ratios, neg_ratios + + +def update_lambdas(const double[:, :, :, ::1] Y, + double[:, :, ::1] X, + double[:, :, :, ::1] X_sigma, + np.ndarray[double, ndim=2, mode='c'] lmbda, + np.ndarray[double, ndim=3, mode='c'] lmbda_sigma, + double[:] lmbda0_eta, + np.ndarray[double, ndim=2, mode='c'] lmbda_eta1, + np.ndarray[double, ndim=3, mode='c'] lmbda_eta2, + double[:, ::1] delta, + double[:, ::1] delta_sigma, + double lmbda_var_prior, + double lmbda_logit_prior, + int max_edges, + int n_neg_subsamples, + double rho, + object rng): + cdef size_t k, p + cdef size_t n_layers = Y.shape[0] + cdef size_t n_features = lmbda.shape[1] + + cdef np.ndarray[np.int64_t, ndim=3, mode='c'] pos_subsamples + cdef np.ndarray[np.int64_t, ndim=3, mode='c'] neg_subsamples + cdef np.ndarray[double, ndim=2, mode='c'] pos_ratios + cdef np.ndarray[double, ndim=2, mode='c'] neg_ratios + cdef np.ndarray[double, ndim=3, mode='c'] omega_pos + cdef np.ndarray[double, ndim=3, mode='c'] omega_neg + + cdef double eta, proba + cdef np.ndarray[double, ndim=1, mode='c'] lmbda0_eta_new = np.zeros( + n_features, dtype=np.float64) + cdef np.ndarray[double, ndim=2, mode='c'] lmbda_new = np.zeros( + (n_layers, n_features), dtype=np.float64) + cdef np.ndarray[double, ndim=3, mode='c'] lmbda_sigma_new = np.zeros( + (n_layers, n_features, n_features), dtype=np.float64) + cdef np.ndarray[double, ndim=1, mode='c'] eta1 + cdef np.ndarray[double, ndim=2, mode='c'] eta2 + cdef np.ndarray[double, ndim=2, mode='c'] lmbda_eta1_new = np.zeros( + (n_layers, n_features), dtype=np.float64) + cdef np.ndarray[double, ndim=3, mode='c'] lmbda_eta2_new = np.zeros( + (n_layers, n_features, n_features), dtype=np.float64) + + # sample dyads + # XXX: this may be slow? + pos_subsamples, neg_subsamples, pos_ratios, neg_ratios = ( + subsample_dyads(np.asarray(Y), max_edges, n_neg_subsamples, rng)) + + # sample local variables + omega_pos, omega_neg = sample_omegas( + X, X_sigma, lmbda, lmbda_sigma, delta, delta_sigma, + pos_subsamples, neg_subsamples) + + # start by updating signs of the reference layer + for p in range(n_features): + eta = calculate_natural_parameters_reference( + Y, X, X_sigma, delta, lmbda[0], omega_pos, omega_neg, + pos_subsamples, neg_subsamples, pos_ratios, neg_ratios, p) + + lmbda0_eta_new[p] = (1 - rho) * lmbda0_eta[p] + rho * eta + + proba = expit(lmbda0_eta_new[p] + lmbda_logit_prior) + lmbda_new[0, p] = 2 * proba - 1 + lmbda_sigma_new[0, p, p] = 1 - lmbda[0, p] ** 2 + + for k in range(1, n_layers): + eta1, eta2 = calculate_natural_parameters( + Y, X, X_sigma, delta, lmbda_var_prior, omega_pos, omega_neg, + pos_subsamples, neg_subsamples, pos_ratios, neg_ratios, k) + + lmbda_eta1_new[k] = (1 - rho) * lmbda_eta1[k] + rho * eta1 + lmbda_eta2_new[k] = (1 - rho) * lmbda_eta2[k] + rho * eta2 + + lmbda_sigma_new[k] = np.linalg.pinv(lmbda_eta2_new[k]) + lmbda_new[k] = np.dot(lmbda_sigma_new[k], lmbda_eta1_new[k]) + + return LambdaParams(lmbda=lmbda_new, lmbda_sigma=lmbda_sigma_new, + lmbda0_eta=lmbda0_eta_new, + lmbda_eta1=lmbda_eta1_new, lmbda_eta2=lmbda_eta2_new) diff --git a/multidynet/node_svi.pyx b/multidynet/node_svi.pyx new file mode 100644 index 0000000..038ea7e --- /dev/null +++ b/multidynet/node_svi.pyx @@ -0,0 +1,470 @@ +# encoding: utf-8 +# cython: language_level=3 +# cython: cdivision=True +# cython: boundscheck=False +# cython: wraparound=False +# cython: nonecheck=False +# cython: initializedcheck=False + +import numpy as np +cimport numpy as np + +from libc.math cimport pow + +from .omega cimport update_omega_svi + + +ctypedef np.int64_t INT + + +class NodeParams(object): + def __init__(self, X, X_sigma, X_cross_cov, X_eta1, X_eta2, + delta, delta_sigma, delta_eta1, delta_eta2): + self.X = X + self.X_sigma = X_sigma + self.X_cross_cov = X_cross_cov + self.X_eta1 = X_eta1 + self.X_eta2 = X_eta2 + self.delta = delta + self.delta_sigma = delta_sigma + self.delta_eta1 = delta_eta1 + self.delta_eta2 = delta_eta2 + + +def sample_omegas(double[:, :, ::1] X, + double[:, :, :, ::1] X_sigma, + double[:, ::1] lmbda, + double[:, :, ::1] lmbda_sigma, + double[:, ::1] delta, + double[:, ::1] delta_sigma, + INT[:, :, ::1] pos_subsamples, + INT[:, :, ::1] neg_subsamples, + int i): + cdef size_t t, k, l + cdef int j + cdef size_t n_time_steps = X.shape[0] + cdef size_t n_layers = lmbda.shape[0] + cdef size_t n_pos_subsamples = pos_subsamples.shape[2] + cdef size_t n_neg_subsamples = neg_subsamples.shape[2] + + cdef np.ndarray[double, ndim=3, mode='c'] omega_pos = np.zeros( + (n_layers, n_time_steps, n_pos_subsamples), np.float64) + cdef np.ndarray[double, ndim=3, mode='c'] omega_neg = np.zeros( + (n_layers, n_time_steps, n_neg_subsamples), np.float64) + + for k in range(n_layers): + for t in range(n_time_steps): + + # loop over positive dyads + for l in range(n_pos_subsamples): + j = pos_subsamples[k, t, l] + if j < 0: + break + + omega_pos[k, t, l] = update_omega_svi( + X[t, i], X_sigma[t, i], X[t, j], X_sigma[t, j], + lmbda[k], lmbda_sigma[k], delta[k, i], delta[k, j], + delta_sigma[k, i], delta_sigma[k, j]) + + # loop over positive dyads + for l in range(n_neg_subsamples): + j = neg_subsamples[k, t, l] + if j < 0: + break + + omega_neg[k, t, l] = update_omega_svi( + X[t, i], X_sigma[t, i], X[t, j], X_sigma[t, j], + lmbda[k], lmbda_sigma[k], delta[k, i], delta[k, j], + delta_sigma[k, i], delta_sigma[k, j]) + + return omega_pos, omega_neg + + +def calculate_natural_parameters_lds(const double[:, :, :, ::1] Y, + double[:, :, ::1] X, + double[:, :, :, ::1] X_sigma, + double[:, ::1] lmbda, + double[:, :, ::1] lmbda_sigma, + double[:, ::1] delta, + double[:, :, ::1] omega_pos, + double[:, :, ::1] omega_neg, + INT[:, :, ::1] pos_subsamples, + INT[:, :, ::1] neg_subsamples, + double[:, ::1] pos_ratios, + double[:, ::1] neg_ratios, + int i): + cdef size_t t, k, l, p, q + cdef int j + cdef size_t n_layers = Y.shape[0] + cdef size_t n_time_steps = Y.shape[1] + cdef size_t n_nodes = Y.shape[2] + cdef size_t n_features = lmbda.shape[1] + cdef size_t n_pos_subsamples = pos_subsamples.shape[2] + cdef size_t n_neg_subsamples = neg_subsamples.shape[2] + cdef double ratio + + cdef double[:, ::1] eta1 = np.zeros( + (n_time_steps, n_features), dtype=np.float64) + cdef double[:, :, ::1] eta2 = np.zeros( + (n_time_steps, n_features, n_features), dtype=np.float64) + + # sample nodes + for t in range(n_time_steps): + for k in range(n_layers): + # loop over positive dyads + ratio = pos_ratios[k, t] + for l in range(n_pos_subsamples): + # extract subample and calculate omega + j = pos_subsamples[k, t, l] + if j < 0: + break + + for p in range(n_features): + eta1[t, p] += ratio * ( + lmbda[k, p] * X[t, j, p] * ( + Y[k, t, i, j] - 0.5 - + omega_pos[k, t, l] * (delta[k, i] + delta[k, j]))) + + for q in range(p + 1): + eta2[t, p, q] += ratio * omega_pos[k, t, l] * ( + (lmbda_sigma[k, p, q] + + lmbda[k, p] * lmbda[k, q]) * + (X_sigma[t, j, p, q] + + X[t, j, p] * X[t, j, q])) + eta2[t, q, p] = eta2[t, p, q] + + # loop over negative (i.e., zero) dyads + ratio = neg_ratios[k, t] + for l in range(n_neg_subsamples): + # extract subample and calculate omega + j = neg_subsamples[k, t, l] + if j < 0: + break + + for p in range(n_features): + eta1[t, p] += ratio * ( + lmbda[k, p] * X[t, j, p] * ( + Y[k, t, i, j] - 0.5 - + omega_neg[k, t, l] * (delta[k, i] + delta[k, j]))) + + for q in range(p + 1): + eta2[t, p, q] += ratio * omega_neg[k, t, l] * ( + (lmbda_sigma[k, p, q] + + lmbda[k, p] * lmbda[k, q]) * + (X_sigma[t, j, p, q] + + X[t, j, p] * X[t, j, q])) + eta2[t, q, p] = eta2[t, p, q] + + return np.asarray(eta1), np.asarray(eta2) + + +def kalman_filter(np.ndarray[double, ndim=2, mode='c'] A, + np.ndarray[double, ndim=3, mode='c'] B, + double tau_prec, + double sigma_prec): + cdef size_t t + cdef size_t n_time_steps = A.shape[0] + cdef size_t n_features = A.shape[1] + + # allocate temporary arrays + cdef np.ndarray[double, ndim=2, mode='c'] mu = np.zeros( + (n_time_steps, n_features)) + cdef np.ndarray[double, ndim=3, mode='c'] sigma = np.zeros( + (n_time_steps, n_features, n_features)) + cdef np.ndarray[double, ndim=3, mode='c'] sigma_inv = np.zeros( + (n_time_steps, n_features, n_features)) + cdef np.ndarray[double, ndim=3, mode='c'] sigma_star = np.zeros( + (n_time_steps - 1, n_features, n_features)) + cdef np.ndarray[double, ndim=2, mode='c'] F_init_inv = ( + tau_prec * np.eye(n_features)) + cdef np.ndarray[double, ndim=2, mode='c'] F_inv = ( + sigma_prec * np.eye(n_features)) + + # t = 1 + sigma_inv[0] = F_init_inv + B[0] + sigma[0] = np.linalg.pinv(sigma_inv[0]) + mu[0] = np.dot(sigma[0], A[0]) + + for t in range(1, n_time_steps): + sigma_star[t-1] = np.linalg.pinv(F_inv + sigma_inv[t-1]) + sigma_inv[t] = F_inv + B[t] - (sigma_prec ** 2) * sigma_star[t-1] + sigma[t] = np.linalg.pinv(sigma_inv[t]) + mu[t] = np.dot(sigma[t], A[t] + + sigma_prec * np.dot( + sigma_star[t-1], np.dot(sigma_inv[t-1], mu[t-1]))) + + return mu, sigma, sigma_inv, sigma_star + + +def kalman_smoother(np.ndarray[double, ndim=2, mode='c'] A, + np.ndarray[double, ndim=3, mode='c'] B, + double tau_prec, + double sigma_prec): + cdef size_t t + cdef size_t n_time_steps = A.shape[0] + cdef size_t n_features = A.shape[1] + + # Allocate temporary arrays + cdef np.ndarray[double, ndim=2, mode='c'] mu + cdef np.ndarray[double, ndim=3, mode='c'] sigma + cdef np.ndarray[double, ndim=3, mode='c'] sigma_inv + cdef np.ndarray[double, ndim=3, mode='c'] sigma_star + cdef np.ndarray[double, ndim=2, mode='c'] eta = np.zeros( + (n_time_steps, n_features)) + cdef np.ndarray[double, ndim=3, mode='c'] psi = np.zeros( + (n_time_steps, n_features, n_features)) + cdef np.ndarray[double, ndim=3, mode='c'] psi_inv = np.zeros( + (n_time_steps, n_features, n_features)) + cdef np.ndarray[double, ndim=3, mode='c'] psi_star = np.zeros( + (n_time_steps, n_features, n_features)) + cdef np.ndarray[double, ndim=2, mode='c'] F_init_inv = ( + tau_prec * np.eye(n_features)) + cdef np.ndarray[double, ndim=2, mode='c'] F_inv = ( + sigma_prec * np.eye(n_features)) + cdef np.ndarray[double, ndim=2, mode='c'] mean = np.zeros( + (n_time_steps, n_features)) + cdef np.ndarray[double, ndim=3, mode='c'] cov = np.zeros( + (n_time_steps, n_features, n_features)) + cdef np.ndarray[double, ndim=3, mode='c'] cross_cov = np.zeros( + (n_time_steps - 1, n_features, n_features)) + + # run the filter for the forward message variables + mu, sigma, sigma_inv, sigma_star = kalman_filter(A, B, tau_prec, sigma_prec) + + # run the smoother + mean[n_time_steps - 1] = mu[n_time_steps - 1] + cov[n_time_steps - 1] = sigma[n_time_steps - 1] + for t in range(n_time_steps - 1, 0, -1): + psi_star[t] = np.linalg.pinv(F_inv + B[t] + psi_inv[t]) + psi_inv[t-1] = F_inv - (sigma_prec ** 2) * psi_star[t] + psi[t-1] = np.linalg.pinv(psi_inv[t-1]) + eta[t-1] = sigma_prec * np.dot(psi[t-1], np.dot(psi_star[t], + A[t] + np.dot(psi_inv[t], eta[t]))) + + # update marginals and cross-covariances + cov[t-1] = np.linalg.pinv(sigma_inv[t-1] + psi_inv[t-1]) + mean[t-1] = np.dot(cov[t-1], + np.dot(sigma_inv[t-1], mu[t-1]) + + np.dot(psi_inv[t-1], eta[t-1])) + cross_cov[t-1] = sigma_prec * np.dot(sigma_star[t-1], + np.linalg.pinv( + F_inv + B[t] + psi_inv[t] - # psi_star^{-1} + (sigma_prec ** 2) * sigma_star[t-1])) + + return mean, cov, cross_cov + + +def calculate_natural_parameters_delta(const double[:, :, :, ::1] Y, + double[:, :, ::1] X, + double[:, ::1] lmbda, + double[:, ::1] delta, + double[:, :, ::1] omega_pos, + double[:, :, ::1] omega_neg, + double delta_var_prior, + INT[:, :, ::1] pos_subsamples, + INT[:, :, ::1] neg_subsamples, + double[:, ::1] pos_ratios, + double[:, ::1] neg_ratios, + int k, + int i): + cdef size_t t, p + cdef int j + cdef size_t n_time_steps = X.shape[0] + cdef size_t n_nodes = X.shape[1] + cdef size_t n_features = X.shape[2] + cdef size_t n_pos_subsamples = pos_subsamples.shape[2] + cdef size_t n_neg_subsamples = neg_subsamples.shape[2] + + cdef double tmp = 0. + cdef double eta1 = 0. + cdef double eta2 = 1. / delta_var_prior + cdef double ratio + + for t in range(n_time_steps): + # positive dyads + for l in range(n_pos_subsamples): + j = pos_subsamples[k, t, l] + if j < 0: + break + + ratio = pos_ratios[k, t] + + eta2 += ratio * omega_pos[k, t, l] + + tmp = delta[k, j] + for p in range(n_features): + tmp += lmbda[k, p] * X[t, i, p] * X[t, j, p] + eta1 += ratio * (Y[k, t, i, j] - 0.5 - omega_pos[k, t, l] * tmp) + + # negative dyads + for l in range(n_neg_subsamples): + j = neg_subsamples[k, t, l] + if j < 0: + break + + ratio = neg_ratios[k, t] + + eta2 += ratio * omega_pos[k, t, l] + + tmp = delta[k, j] + for p in range(n_features): + tmp += lmbda[k, p] * X[t, i, p] * X[t, j, p] + eta1 += ratio * (Y[k, t, i, j] - 0.5 - omega_neg[k, t, l] * tmp) + + return eta1, eta2 + + +def subsample_nodes(np.ndarray[double, ndim=4, mode='c'] Y, + int max_degree, + int n_neg_multiplier, + int i, + object rng): + cdef size_t n_layers = Y.shape[0] + cdef size_t n_time_steps = Y.shape[1] + cdef size_t n_nodes = Y.shape[2] + cdef size_t k, t + cdef size_t n_pos_samples, n_neg_samples, n_neg_subsamples + + cdef np.ndarray[double, ndim=1, mode='c'] y_row + cdef np.ndarray[INT, ndim=1, mode='c'] indices + cdef np.ndarray[INT, ndim=3, mode='c'] pos_subsamples = np.full( + (n_layers, n_time_steps, max_degree), -1, dtype=np.int64) + cdef np.ndarray[double, ndim=2, mode='c'] pos_ratios = np.empty( + (n_layers, n_time_steps), dtype=np.float64) + cdef np.ndarray[INT, ndim=3, mode='c'] neg_subsamples = np.full( + (n_layers, n_time_steps, n_neg_multiplier * max_degree), -1, dtype=np.int64) + cdef np.ndarray[double, ndim=2, mode='c'] neg_ratios = np.empty( + (n_layers, n_time_steps), dtype=np.float64) + + for k in range(n_layers): + for t in range(n_time_steps): + y_row = Y[k, t, i].copy() + + # XXX: positive / negative dyads could be cached + # include all edges (positive dyads) + indices = np.where(y_row == 1)[0] + n_pos_samples = indices.shape[0] + if n_pos_samples > 0: + #indices = rng.choice(indices, size=n_pos_subsamples) + pos_subsamples[k, t, :n_pos_samples] = indices + pos_ratios[k, t] = n_pos_samples / ( indices.shape[0]) + + # XXX: postive / negative dyads could be cached + # sample zero dyads (negative connections) + y_row[i] = 1 # so we do not sample the diagonal + indices = np.where(y_row == 0)[0] + n_neg_samples = indices.shape[0] + if n_neg_samples: + n_neg_subsamples = min(n_neg_samples, n_neg_multiplier * n_pos_samples) + if n_neg_subsamples < n_neg_samples: + indices = rng.choice(indices, size=n_neg_subsamples) + neg_subsamples[k, t, :indices.shape[0]] = indices + neg_ratios[k, t] = n_neg_samples / ( indices.shape[0]) + + return pos_subsamples, neg_subsamples, pos_ratios, neg_ratios + + +def update_node_effects(const double[:, :, :, ::1] Y, + np.ndarray[double, ndim=3, mode='c'] X, + np.ndarray[double, ndim=4, mode='c'] X_sigma, + np.ndarray[double, ndim=3, mode='c'] X_eta1, + np.ndarray[double, ndim=4, mode='c'] X_eta2, + double[:, ::1] lmbda, + double[:, :, ::1] lmbda_sigma, + double[:, ::1] delta, + double[:, ::1] delta_sigma, + double[:, ::1] delta_eta1, + double[:, ::1] delta_eta2, + double delta_var_prior, + double tau_prec, + double sigma_prec, + int max_degree, + int n_neg_multiplier, + double rho, + object rng): + cdef size_t i, k + cdef size_t n_layers = Y.shape[0] + cdef size_t n_time_steps = X.shape[0] + cdef size_t n_nodes = X.shape[1] + cdef size_t n_features = X.shape[2] + + # local variables + cdef np.ndarray[INT, ndim=3, mode='c'] pos_subsamples + cdef np.ndarray[INT, ndim=3, mode='c'] neg_subsamples + cdef np.ndarray[double, ndim=2, mode='c'] pos_ratios + cdef np.ndarray[double, ndim=2, mode='c'] neg_ratios + cdef np.ndarray[double, ndim=3, mode='c'] omega_pos + cdef np.ndarray[double, ndim=3, mode='c'] omega_neg + + # updated natural parameters (lds) + cdef np.ndarray[double, ndim=2, mode='c'] lds_eta1 + cdef np.ndarray[double, ndim=3, mode='c'] lds_eta2 + cdef np.ndarray[double, ndim=3, mode='c'] X_eta1_new = np.zeros( + (n_nodes, n_time_steps, n_features), dtype=np.float64) + cdef np.ndarray[double, ndim=4, mode='c'] X_eta2_new = np.zeros( + (n_nodes, n_time_steps, n_features, n_features), dtype=np.float64) + + # updated natural parameters (sociability) + cdef double eta1 + cdef double eta2 + cdef double[:, ::1] delta_eta1_new = np.zeros( + (n_layers, n_nodes), dtype=np.float64) + cdef double[:, ::1] delta_eta2_new = np.zeros( + (n_layers, n_nodes), dtype=np.float64) + + # updated moments (lds) + cdef np.ndarray[double, ndim=3, mode='c'] X_new = np.zeros( + (n_time_steps, n_nodes, n_features), dtype=np.float64) + cdef np.ndarray[double, ndim=4, mode='c'] X_sigma_new = np.zeros( + (n_time_steps, n_nodes, n_features, n_features), dtype=np.float64) + cdef np.ndarray[double, ndim=4, mode='c'] X_cross_cov_new = np.zeros( + (n_time_steps - 1, n_nodes, n_features, n_features), dtype=np.float64) + + # updated moments (sociability) + cdef double[:, ::1] delta_new = np.zeros( + (n_layers, n_nodes), dtype=np.float64) + cdef double[:, ::1] delta_sigma_new = np.zeros( + (n_layers, n_nodes), dtype=np.float64) + + for i in range(n_nodes): + # XXX: this may be slow? + # determine subsamples + pos_subsamples, neg_subsamples, pos_ratios, neg_ratios = ( + subsample_nodes( + np.asarray(Y), max_degree, n_neg_multiplier, i, rng)) + + # sample local variables + omega_pos, omega_neg = sample_omegas( + X, X_sigma, lmbda, lmbda_sigma, delta, delta_sigma, + pos_subsamples, neg_subsamples, i) + + # calculate natural parameters on this subsample + lds_eta1, lds_eta2 = calculate_natural_parameters_lds( + Y, X, X_sigma, lmbda, lmbda_sigma, delta, omega_pos, omega_neg, + pos_subsamples, neg_subsamples, pos_ratios, neg_ratios, i) + + # take a step along the gradient + X_eta1_new[i] = (1 - rho) * X_eta1[i] + rho * lds_eta1 + X_eta2_new[i] = (1 - rho) * X_eta2[i] + rho * lds_eta2 + + # update latent position expectations + X_new[:, i], X_sigma_new[:, i], X_cross_cov_new[:, i] = kalman_smoother( + X_eta1_new[i], X_eta2_new[i], tau_prec, sigma_prec) + + # update sociability parameters + for k in range(n_layers): + eta1, eta2 = calculate_natural_parameters_delta( + Y, X, lmbda, delta, omega_pos, omega_neg, + delta_var_prior, pos_subsamples, neg_subsamples, + pos_ratios, neg_ratios, k, i) + + delta_eta1_new[k, i] = (1 - rho) * delta_eta1[k, i] + rho * eta1 + delta_eta2_new[k, i] = (1 - rho) * delta_eta2[k, i] + rho * eta2 + + delta_sigma_new[k, i] = 1. / delta_eta2_new[k, i] + delta_new[k, i] = delta_sigma_new[k, i] * delta_eta1_new[k, i] + + return NodeParams(X=X_new, X_sigma=X_sigma_new, X_cross_cov=X_cross_cov_new, + X_eta1=X_eta1_new, X_eta2=X_eta2_new, + delta=delta_new, delta_sigma=delta_sigma_new, + delta_eta1=delta_eta1_new, delta_eta2=delta_eta2_new) diff --git a/multidynet/omega.pxd b/multidynet/omega.pxd new file mode 100644 index 0000000..1c0202e --- /dev/null +++ b/multidynet/omega.pxd @@ -0,0 +1,13 @@ +# cython: language_level=3 + +# update a single auxilary variable +cdef double update_omega_svi(double[::1] Xit, + double[:, ::1] Xit_sigma, + double[::1] Xjt, + double[:, ::1] Xjt_sigma, + double[::1] lmbdak, + double[:, ::1] lmbdak_sigma, + double deltaki, + double deltakj, + double deltaki_sigma, + double deltakj_sigma) diff --git a/multidynet/omega.pyx b/multidynet/omega.pyx index b5ead8a..c3b5edc 100644 --- a/multidynet/omega.pyx +++ b/multidynet/omega.pyx @@ -12,18 +12,55 @@ import numpy as np cimport numpy as np -def update_omega_single(double interceptk, - double interceptk_sigma, - double[::1] Xit, - double[:, ::1] Xit_sigma, - double[::1] Xjt, - double[:, ::1] Xjt_sigma, - double[::1] lmbdak, - double[:, ::1] lmbdak_sigma, - double deltaki, - double deltakj, - double deltaki_sigma, - double deltakj_sigma): +cdef double update_omega_svi(double[::1] Xit, + double[:, ::1] Xit_sigma, + double[::1] Xjt, + double[:, ::1] Xjt_sigma, + double[::1] lmbdak, + double[:, ::1] lmbdak_sigma, + double deltaki, + double deltakj, + double deltaki_sigma, + double deltakj_sigma): + cdef double psi_sq = 0. + cdef double c_omega = 0. + cdef double omega = 0. + cdef size_t n_features = Xit.shape[0] + cdef int p, q = 0 + + # calculate the natural parameter + psi_sq += deltaki_sigma + deltaki ** 2 + psi_sq += deltakj_sigma + deltakj ** 2 + psi_sq += 2 * deltaki * deltakj + for p in range(n_features): + psi_sq += (2 * (deltaki + deltakj) * + lmbdak[p] * Xit[p] * Xjt[p]) + + for q in range(n_features): + psi_sq += ((lmbdak_sigma[p, q] + lmbdak[p] * lmbdak[q]) * + (Xit_sigma[p, q] + Xit[p] * Xit[q]) * + (Xjt_sigma[p, q] + Xjt[p] * Xjt[q])) + + # calculate mean of a PG(1, sqrt(c_omega)) random variable + c_omega = sqrt(psi_sq) + omega = tanh(0.5 * c_omega) + omega /= (2. * c_omega) + + return omega + + +cdef update_omega_single(double interceptk, + double interceptk_sigma, + double[::1] Xit, + double[:, ::1] Xit_sigma, + double[::1] Xjt, + double[:, ::1] Xjt_sigma, + double[::1] lmbdak, + double[:, ::1] lmbdak_sigma, + double deltaki, + double deltakj, + double deltaki_sigma, + double deltakj_sigma): cdef double psi_sq = 0. cdef double c_omega = 0. cdef double omega = 0. diff --git a/multidynet/svi.py b/multidynet/svi.py new file mode 100644 index 0000000..95b3114 --- /dev/null +++ b/multidynet/svi.py @@ -0,0 +1,385 @@ +import numbers +import warnings + +import numpy as np +import scipy.sparse as sp + +from joblib import Parallel, delayed +from scipy.special import logit, gammainc, expit +from sklearn.exceptions import ConvergenceWarning +from sklearn.metrics import roc_auc_score +from sklearn.utils import check_array, check_random_state +from sklearn.linear_model import LogisticRegression +from tqdm import tqdm + +from .node_svi import update_node_effects +from .lmbdas_svi import update_lambdas +from .variances import update_tau_sq, update_sigma_sq +from .log_likelihood import log_likelihood +from .metrics import calculate_auc + + +__all__ = ['DynamicMultilayerNetworkLSM'] + + + +class ModelParameters(object): + def __init__(self, X, X_sigma, X_cross_cov, X_eta1, + X_eta2, lmbda, lmbda_sigma, + lmbda_logit_prior, lmbda0_eta, lmbda_eta1, lmbda_eta2, + delta, delta_sigma, delta_eta1, delta_eta2, + a_tau_sq, b_tau_sq, c_sigma_sq, + d_sigma_sq, logp=None): + self.X_ = X + self.X_sigma_ = X_sigma + self.X_cross_cov_ = X_cross_cov + self.X_eta1_ = X_eta1 + self.X_eta2_ = X_eta2 + + self.lambda_ = lmbda + self.lambda_sigma_ = lmbda_sigma + self.lambda_logit_prior_ = lmbda_logit_prior + self.lambda0_eta_ = lmbda0_eta + self.lambda_eta1_ = lmbda_eta1 + self.lambda_eta2_ = lmbda_eta2 + + self.delta_ = delta + self.delta_sigma_ = delta_sigma + self.delta_eta1_ = delta_eta1 + self.delta_eta2_ = delta_eta2 + + self.a_tau_sq_ = a_tau_sq + self.b_tau_sq_ = b_tau_sq + self.c_sigma_sq_ = c_sigma_sq + self.d_sigma_sq_ = d_sigma_sq + self.converged_ = False + + self.logp_ = [] if logp is None else logp + + +def initialize_node_effects_single(Y): + n_time_steps, n_nodes, _ = Y.shape + + n_dyads = int(0.5 * n_nodes * (n_nodes - 1)) + dyads = np.tril_indices_from(Y[0], k=-1) + y_vec = np.zeros(n_time_steps * n_dyads) + + # construct dummy node indicators + cols = np.r_[dyads[0], dyads[1]] + rows = np.r_[np.arange(n_dyads), np.arange(n_dyads)] + x_dummy = sp.coo_matrix((np.ones(2 * n_dyads), (rows, cols)), + shape=(n_dyads, n_nodes)) + + # dyad target + for t in range(n_time_steps): + yt_vec = Y[t][dyads] + y_vec[(t * n_dyads):((t+1) * n_dyads)] = Y[t][dyads] + if t > 0: + X = sp.vstack((X, x_dummy)) + else: + X = x_dummy.copy() + X = X.tocsr() + + # remove missing values + non_missing = y_vec != -1.0 + + logreg = LogisticRegression(fit_intercept=False, C=1e5) + logreg.fit(X[non_missing], y_vec[non_missing]) + + return logreg.coef_[0] + + +def initialize_node_effects(Y): + n_layers, n_time_steps, n_nodes, _ = Y.shape + + delta = np.zeros((n_layers, n_nodes)) + for k in range(n_layers): + delta[k] = initialize_node_effects_single(Y[k]) + + return delta + + +def initialize_parameters(Y, n_features, lambda_odds_prior, lambda_var_prior, + delta_var_prior, + a, b, c, d, random_state): + rng = check_random_state(random_state) + + n_layers, n_time_steps, n_nodes, _ = Y.shape + + # intialize latent space randomly + X = rng.randn(n_time_steps, n_nodes, n_features) + + # intialize to marginal covariances + sigma_init = np.eye(n_features) + X_sigma = np.tile( + sigma_init[None, None], reps=(n_time_steps, n_nodes, 1, 1)) + + # initialize cross-covariances + cross_init = np.eye(n_features) + X_cross_cov = np.tile( + cross_init[None, None], reps=(n_time_steps - 1, n_nodes, 1, 1)) + + # initialize natural parameters + X_eta1 = np.zeros((n_nodes, n_time_steps, n_features), dtype=np.float64) + X_eta2 = np.zeros( + (n_nodes, n_time_steps, n_features, n_features), dtype=np.float64) + + # initialize node-effects based on a logistic regression fit + delta = initialize_node_effects(Y) + delta_sigma = delta_var_prior * np.ones((n_layers, n_nodes)) + + # initialize natural parameters + delta_eta1 = np.zeros((n_layers, n_nodes), dtype=np.float64) + delta_eta2 = np.zeros((n_layers, n_nodes), dtype=np.float64) + + # intialize to prior means + lmbda = np.sqrt(2) * rng.randn(n_layers, n_features) + lmbda[0] = ( + 2 * (lambda_odds_prior / (1. + lambda_odds_prior)) - 1) + lmbda_sigma = lambda_var_prior * np.ones( + (n_layers, n_features, n_features)) + lmbda_sigma[0] = ( + (1 - lmbda[0, 0] ** 2) * np.eye(n_features)) + lmbda_logit_prior = np.log(lambda_odds_prior) + + # initialize natural parameters + lmbda0_eta = np.zeros(n_features, dtype=np.float64) + lmbda_eta1 = np.zeros((n_layers, n_features), dtype=np.float64) + lmbda_eta2 = np.zeros((n_layers, n_features, n_features), dtype=np.float64) + + # initialize based on prior information + a_tau_sq = a + b_tau_sq = b + c_sigma_sq = c + d_sigma_sq = d + + return ModelParameters( + X=X, X_sigma=X_sigma, X_cross_cov=X_cross_cov, + X_eta1=X_eta1, X_eta2=X_eta2, + lmbda=lmbda, + lmbda_sigma=lmbda_sigma, lmbda_logit_prior=lmbda_logit_prior, + lmbda0_eta=lmbda0_eta, lmbda_eta1=lmbda_eta1, lmbda_eta2=lmbda_eta2, + delta=delta, delta_sigma=delta_sigma, delta_eta1=delta_eta1, + delta_eta2=delta_eta2, + a_tau_sq=a_tau_sq, b_tau_sq=b_tau_sq, c_sigma_sq=c_sigma_sq, + d_sigma_sq=d_sigma_sq) + + + +def optimize_elbo(Y, n_features, lambda_odds_prior, lambda_var_prior, + delta_var_prior, tau_sq, sigma_sq, a, b, c, d, + gamma, max_degree, max_edges, n_neg_multiplier, + max_iter, tol, random_state, verbose=True): + + # initialize parameters of the model + random_state = check_random_state(random_state) + + model = initialize_parameters( + Y, n_features, lambda_odds_prior, lambda_var_prior, delta_var_prior, + a, b, c, d, random_state) + + for n_iter in tqdm(range(max_iter), disable=not verbose): + # stochastic gradient descent + + # update step size + rho = (1 + n_iter) ** (-gamma) + + # latent trajectory updates + tau_sq_prec = ( + model.a_tau_sq_ / model.b_tau_sq_ if tau_sq == 'auto' else + 1. / tau_sq) + sigma_sq_prec = ( + model.c_sigma_sq_ / model.d_sigma_sq_ if sigma_sq == 'auto' else + 1. / sigma_sq) + + # update latent positions and degree effects + node_params = update_node_effects( + Y, model.X_, model.X_sigma_, model.X_eta1_, model.X_eta2_, + model.lambda_, model.lambda_sigma_, model.delta_, + model.delta_sigma_, model.delta_eta1_, model.delta_eta2_, + delta_var_prior, tau_sq_prec, sigma_sq_prec, + max_degree, n_neg_multiplier, rho, random_state) + + # update lambda values + lmbda_params = update_lambdas( + Y, model.X_, model.X_sigma_, model.lambda_, + model.lambda_sigma_, model.lambda0_eta_, model.lambda_eta1_, + model.lambda_eta2_, model.delta_, model.delta_sigma_, + lambda_var_prior, model.lambda_logit_prior_, + max_edges, n_neg_multiplier, rho, random_state) + + # update intial variance of the latent space + if tau_sq == 'auto': + a_tau_sq, b_tau_sq = update_tau_sq( + Y, model.X_, model.X_sigma_, a, b) + + a_tau_sq = (1 - rho) * model.a_tau_sq_ + rho * a_tau_sq + b_tau_sq = (1 - rho) * model.b_tau_sq_ + rho * b_tau_sq + else: + a_tau_sq = model.a_tau_sq_ + b_tau_sq = model.b_tau_sq_ + + # update step sizes + if sigma_sq == 'auto': + c_sigma_sq, d_sigma_sq = update_sigma_sq( + Y, model.X_, model.X_sigma_, model.X_cross_cov_, c, d) + + c_sigma_sq = (1 - rho) * model.c_sigma_sq_ + rho * c_sigma_sq + d_sigma_sq = (1 - rho) * model.d_sigma_sq_ + rho * d_sigma_sq + else: + c_sigma_sq = model.c_sigma_sq_ + d_sigma_sq = model.d_sigma_sq_ + + + # XXX: calculate change in latent space parameters + change = ((model.X_ - node_params.X) ** 2).sum() + model.logp_.append(change) + + # store updated model parameters + model = ModelParameters( + X=node_params.X, X_sigma=node_params.X_sigma, + X_cross_cov=node_params.X_cross_cov, X_eta1=node_params.X_eta1, + X_eta2=node_params.X_eta2, lmbda=lmbda_params.lmbda, + lmbda_logit_prior=model.lambda_logit_prior_, + lmbda_sigma=lmbda_params.lmbda_sigma, + lmbda0_eta=lmbda_params.lmbda0_eta, + lmbda_eta1=lmbda_params.lmbda_eta1, + lmbda_eta2=lmbda_params.lmbda_eta2, + delta=node_params.delta, + delta_sigma=node_params.delta_sigma, + delta_eta1=node_params.delta_eta1, + delta_eta2=node_params.delta_eta2, a_tau_sq=a_tau_sq, + b_tau_sq=b_tau_sq, c_sigma_sq=c_sigma_sq, d_sigma_sq=d_sigma_sq, + logp=model.logp_) + + # We could calculate the expected log-likelihood on dyads sampled + # to update lmbda + + # check convergence + #change = loglik - prev_loglik + if change < tol: + model.converged_ = True + #model.logp_ = np.asarray(model.logp_) + break + + return model + + +def calculate_probabilities(X, lmbda, delta): + n_layers = lmbda.shape[0] + n_time_steps = X.shape[0] + n_nodes = X.shape[1] + + probas = np.zeros( + (n_layers, n_time_steps, n_nodes, n_nodes), dtype=np.float64) + for k in range(n_layers): + for t in range(n_time_steps): + deltak = delta[k].reshape(-1, 1) + eta = np.add(deltak, deltak.T) + np.dot(X[t] * lmbda[k], X[t].T) + probas[k, t] = expit(eta) + + return probas + + +class DynamicMultilayerNetworkLSM(object): + def __init__(self, n_features=2, + lambda_odds_prior=2, + lambda_var_prior=4, + delta_var_prior=4, + tau_sq='auto', sigma_sq='auto', + a=4.0, b=8.0, c=10, d=0.1, + gamma=0.75, n_neg_multiplier=2, + n_init=1, max_iter=500, tol=1e-2, + n_jobs=-1, random_state=42): + self.n_features = n_features + self.lambda_odds_prior = lambda_odds_prior + self.lambda_var_prior = lambda_var_prior + self.delta_var_prior = delta_var_prior + self.tau_sq = tau_sq + self.sigma_sq = sigma_sq + self.a = a + self.b = b + self.c = c + self.d = d + self.gamma = gamma + self.n_neg_multiplier = n_neg_multiplier + self.n_init = n_init + self.max_iter = max_iter + self.tol = tol + self.n_jobs = n_jobs + self.random_state = random_state + + def fit(self, Y): + """ + Parameters + ---------- + Y : array-like, shape (n_layers, n_time_steps, n_nodes, n_nodes) + """ + Y = check_array(Y, order='C', dtype=np.float64, + ensure_2d=False, allow_nd=True, copy=False) + + n_nodes = Y.shape[2] + + # maximum degree, i.e. the maximum number of edges connected to a node + max_degree = int(np.count_nonzero(Y == 1, axis=3).max()) + + # maximum edge count + max_edges = int((np.count_nonzero(Y == 1, axis=(2, 3)) / 2).max()) + + if self.gamma >= 1.0: + raise ValueError('Decay rate `gamma` must be between 0 and 1') + + random_state = check_random_state(self.random_state) + + # run the elbo optimization over different initializations + seeds = random_state.randint(np.iinfo(np.int32).max, size=self.n_init) + verbose = True if self.n_init == 1 else False + models = Parallel(n_jobs=self.n_jobs)(delayed(optimize_elbo)( + Y, self.n_features, self.lambda_odds_prior, + self.lambda_var_prior, self.delta_var_prior, + self.tau_sq, self.sigma_sq, self.a, self.b, self.c, self.d, + self.gamma, max_degree, max_edges, self.n_neg_multiplier, + self.max_iter, self.tol, seed, verbose=verbose) + for seed in seeds) + + # choose model with the largest convergence criteria + best_model = models[0] + #best_criteria = models[0].logp_[-1] + #for i in range(1, len(models)): + # if models[i].logp_[-1] > best_criteria: + # best_model = models[i] + + if not best_model.converged_: + warnings.warn('Best model did not converge. ' + 'Try a different random initialization, ' + 'or increase max_iter, tol ' + 'or check for degenerate data.', ConvergenceWarning) + + self._set_parameters(best_model) + + # calculate dyad-probabilities + self.probas_ = calculate_probabilities( + self.X_, self.lambda_, self.delta_) + + # calculate in-sample AUC + self.auc_ = calculate_auc(Y, self.probas_) + + return self + + def _set_parameters(self, model): + self.X_ = model.X_ + self.X_sigma_ = model.X_sigma_ + self.X_cross_cov_ = model.X_cross_cov_ + self.delta_ = np.asarray(model.delta_) + self.delta_sigma_ = np.asarray(model.delta_sigma_) + self.lambda_ = model.lambda_ + self.lambda_[0] = np.sign(model.lambda_[0]) + self.lambda_proba_ = (model.lambda_[0] + 1) / 2. + self.lambda_sigma_ = model.lambda_sigma_ + self.a_tau_sq_ = model.a_tau_sq_ + self.b_tau_sq_ = model.b_tau_sq_ + self.tau_sq_ = self.b_tau_sq_ / (self.a_tau_sq_ - 1) + self.c_sigma_sq_ = model.c_sigma_sq_ + self.d_sigma_sq_ = model.d_sigma_sq_ + self.sigma_sq_ = self.d_sigma_sq_ / (self.c_sigma_sq_ - 1) + self.logp_ = np.asarray(model.logp_) diff --git a/notebooks/ICEWS [Obama Administration].ipynb b/notebooks/ICEWS [Obama Administration].ipynb index 1abf3cc..016c31f 100644 --- a/notebooks/ICEWS [Obama Administration].ipynb +++ b/notebooks/ICEWS [Obama Administration].ipynb @@ -22,12 +22,12 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## 4 Layers, 96 time-stamps (monthly), 65 nodes, 20% of nodes held-out of each network" + "## 4 Layers, 96 time-stamps (monthly), 65 nodes, 20% of dyads held-out of each network" ] }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 5, "metadata": {}, "outputs": [], "source": [ @@ -43,7 +43,27 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "579.0" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "(np.count_nonzero(Y_train == 1, axis=(2, 3)) / 2).max()" + ] + }, + { + "cell_type": "code", + "execution_count": 9, "metadata": {}, "outputs": [], "source": [ @@ -8292,6 +8312,606 @@ "\n", "fig.subplots_adjust(right=0.8)" ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Subsamples" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": {}, + "outputs": [], + "source": [ + "from multidynet.batch import subsample_nodes\n", + "\n", + "rng = np.random.RandomState(123)\n", + "\n", + "pos_subsample, neg_subsample, pos_ratio, neg_ratio = subsample_nodes(Y, 3, 20, 10, rng)" + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[0.27027027, 0.22727273, 0.3030303 , 0.25 , 0.25 ,\n", + " 0.29411765, 0.25 , 0.2173913 , 0.21276596, 0.23809524,\n", + " 0.25 , 0.22222222, 0.19607843, 0.24390244, 0.22727273,\n", + " 0.23255814, 0.24390244, 0.2173913 , 0.22222222, 0.22222222,\n", + " 0.25641026, 0.23809524, 0.23255814, 0.22727273, 0.23255814,\n", + " 0.2173913 , 0.2173913 , 0.2173913 , 0.23255814, 0.23809524,\n", + " 0.23255814, 0.19230769, 0.25641026, 0.23255814, 0.24390244,\n", + " 0.19607843, 0.22727273, 0.22727273, 0.2173913 , 0.25641026,\n", + " 0.23809524, 0.25641026, 0.23255814, 0.26315789, 0.22727273,\n", + " 0.22727273, 0.23809524, 0.20408163, 0.2 , 0.2173913 ,\n", + " 0.18867925, 0.21276596, 0.22222222, 0.2173913 , 0.20833333,\n", + " 0.24390244, 0.24390244, 0.22222222, 0.23255814, 0.20408163,\n", + " 0.2 , 0.22727273, 0.25 , 0.2173913 , 0.22727273,\n", + " 0.3125 , 0.23809524, 0.24390244, 0.3030303 , 0.25641026,\n", + " 0.33333333, 0.23255814, 0.23809524, 0.22727273, 0.22222222,\n", + " 0.24390244, 0.24390244, 0.24390244, 0.26315789, 0.23809524,\n", + " 0.25641026, 0.26315789, 0.26315789, 0.2173913 , 0.25 ,\n", + " 0.25 , 0.23809524, 0.22222222, 0.22222222, 0.20408163,\n", + " 0.22727273, 0.21276596, 0.27027027, 0.22727273, 0.25641026,\n", + " 0.2173913 ],\n", + " [0.17241379, 0.16666667, 0.16666667, 0.16129032, 0.16129032,\n", + " 0.16666667, 0.17241379, 0.16949153, 0.18181818, 0.16666667,\n", + " 0.15873016, 0.16666667, 0.16666667, 0.15873016, 0.16129032,\n", + " 0.17241379, 0.18181818, 0.16393443, 0.16129032, 0.16666667,\n", + " 0.16949153, 0.16129032, 0.17857143, 0.16393443, 0.16393443,\n", + " 0.16666667, 0.16129032, 0.16949153, 0.15625 , 0.16129032,\n", + " 0.15625 , 0.16129032, 0.15873016, 0.16129032, 0.15873016,\n", + " 0.17241379, 0.16393443, 0.16393443, 0.16129032, 0.16393443,\n", + " 0.17241379, 0.15625 , 0.16129032, 0.16666667, 0.15873016,\n", + " 0.16949153, 0.15625 , 0.16129032, 0.15625 , 0.17241379,\n", + " 0.15873016, 0.16129032, 0.16393443, 0.15873016, 0.16666667,\n", + " 0.16129032, 0.15873016, 0.15873016, 0.16129032, 0.15625 ,\n", + " 0.16393443, 0.16129032, 0.15625 , 0.15873016, 0.15873016,\n", + " 0.16393443, 0.16129032, 0.16393443, 0.17241379, 0.16949153,\n", + " 0.16949153, 0.16129032, 0.16666667, 0.16129032, 0.16949153,\n", + " 0.16949153, 0.16393443, 0.16393443, 0.16393443, 0.16393443,\n", + " 0.16393443, 0.16393443, 0.16129032, 0.16393443, 0.16393443,\n", + " 0.16949153, 0.16949153, 0.16129032, 0.16393443, 0.16393443,\n", + " 0.16666667, 0.16393443, 0.16129032, 0.16129032, 0.17241379,\n", + " 0.15873016],\n", + " [0.19607843, 0.19230769, 0.17241379, 0.18867925, 0.1754386 ,\n", + " 0.18867925, 0.18867925, 0.16949153, 0.18181818, 0.18181818,\n", + " 0.18181818, 0.1754386 , 0.17857143, 0.19230769, 0.18518519,\n", + " 0.18518519, 0.19607843, 0.18518519, 0.17857143, 0.18867925,\n", + " 0.18181818, 0.17241379, 0.19230769, 0.18867925, 0.18518519,\n", + " 0.17241379, 0.17241379, 0.18867925, 0.1754386 , 0.1754386 ,\n", + " 0.16949153, 0.18867925, 0.17241379, 0.18518519, 0.17857143,\n", + " 0.16129032, 0.16129032, 0.18181818, 0.17241379, 0.1754386 ,\n", + " 0.17857143, 0.18181818, 0.18181818, 0.18181818, 0.18181818,\n", + " 0.18518519, 0.19230769, 0.18518519, 0.16666667, 0.17241379,\n", + " 0.1754386 , 0.18518519, 0.16666667, 0.1754386 , 0.1754386 ,\n", + " 0.1754386 , 0.17241379, 0.18181818, 0.17857143, 0.17857143,\n", + " 0.18181818, 0.18867925, 0.1754386 , 0.1754386 , 0.17857143,\n", + " 0.21276596, 0.19230769, 0.19607843, 0.20833333, 0.18518519,\n", + " 0.19607843, 0.18867925, 0.1754386 , 0.17241379, 0.17857143,\n", + " 0.18867925, 0.18518519, 0.17857143, 0.18867925, 0.19230769,\n", + " 0.18867925, 0.16949153, 0.18518519, 0.17857143, 0.18181818,\n", + " 0.1754386 , 0.18181818, 0.1754386 , 0.1754386 , 0.18181818,\n", + " 0.18518519, 0.17857143, 0.19230769, 0.18867925, 0.18181818,\n", + " 0.17857143],\n", + " [0.19230769, 0.17241379, 0.18181818, 0.16666667, 0.17857143,\n", + " 0.19230769, 0.18181818, 0.17241379, 0.16949153, 0.1754386 ,\n", + " 0.17241379, 0.17857143, 0.1754386 , 0.17241379, 0.17857143,\n", + " 0.17857143, 0.18867925, 0.1754386 , 0.17857143, 0.17241379,\n", + " 0.17857143, 0.16666667, 0.1754386 , 0.18181818, 0.16666667,\n", + " 0.16129032, 0.17241379, 0.1754386 , 0.18518519, 0.16949153,\n", + " 0.16666667, 0.17241379, 0.16666667, 0.16666667, 0.17241379,\n", + " 0.18181818, 0.16393443, 0.17857143, 0.16393443, 0.17241379,\n", + " 0.16949153, 0.17241379, 0.17241379, 0.16393443, 0.16949153,\n", + " 0.16949153, 0.17857143, 0.16949153, 0.1754386 , 0.16949153,\n", + " 0.18867925, 0.17241379, 0.17241379, 0.16666667, 0.16666667,\n", + " 0.1754386 , 0.17241379, 0.16666667, 0.16949153, 0.16666667,\n", + " 0.17857143, 0.18181818, 0.17241379, 0.16666667, 0.18518519,\n", + " 0.17241379, 0.16666667, 0.18518519, 0.1754386 , 0.18518519,\n", + " 0.1754386 , 0.18181818, 0.17857143, 0.19607843, 0.18867925,\n", + " 0.17857143, 0.18181818, 0.17857143, 0.1754386 , 0.18181818,\n", + " 0.18181818, 0.17241379, 0.18518519, 0.17241379, 0.16666667,\n", + " 0.1754386 , 0.1754386 , 0.18181818, 0.1754386 , 0.17857143,\n", + " 0.17241379, 0.1754386 , 0.17241379, 0.17241379, 0.1754386 ,\n", + " 0.1754386 ]])" + ] + }, + "execution_count": 29, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "neg_ratio" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(4, 96, 20)" + ] + }, + "execution_count": 30, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "pos_subsample.shape" + ] + }, + { + "cell_type": "code", + "execution_count": 38, + "metadata": {}, + "outputs": [], + "source": [ + "y = np.zeros((5, 5))\n", + "y[np.tril_indices_from(y, k=-1)] = np.arange(0.5 * 5 * (5- 1)) + 1\n", + "y += y.T\n", + "y = y.astype(np.int)\n", + "y_vec = y[np.tril_indices_from(y, k=-1)]" + ] + }, + { + "cell_type": "code", + "execution_count": 39, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[ 0, 1, 2, 4, 7],\n", + " [ 1, 0, 3, 5, 8],\n", + " [ 2, 3, 0, 6, 9],\n", + " [ 4, 5, 6, 0, 10],\n", + " [ 7, 8, 9, 10, 0]])" + ] + }, + "execution_count": 39, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "y" + ] + }, + { + "cell_type": "code", + "execution_count": 40, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[ 0, 1, 2, 4, 7],\n", + " [ 1, 0, 3, 5, 8],\n", + " [ 2, 3, 0, 6, 9],\n", + " [ 4, 5, 6, 0, 10],\n", + " [ 7, 8, 9, 10, 0]])" + ] + }, + "execution_count": 40, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "y" + ] + }, + { + "cell_type": "code", + "execution_count": 41, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "8" + ] + }, + "execution_count": 41, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "i = 4\n", + "j = 1\n", + "y_vec[int(i * (i - 1) / 2 + j)]" + ] + }, + { + "cell_type": "code", + "execution_count": 66, + "metadata": {}, + "outputs": [], + "source": [ + "from math import ceil, floor\n", + "\n", + "idx = 0\n", + "i = ceil(np.sqrt(2 * (idx + 1) + 0.25) - 0.5)\n", + "j = int(idx - 0.5 * i * (i - 1)) " + ] + }, + { + "cell_type": "code", + "execution_count": 75, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "1 1 (1, 0)\n", + "2 2 (2, 0)\n", + "3 3 (2, 1)\n", + "4 4 (3, 0)\n", + "5 5 (3, 1)\n", + "6 6 (3, 2)\n", + "7 7 (4, 0)\n", + "8 8 (4, 1)\n", + "9 9 (4, 2)\n", + "10 10 (4, 3)\n" + ] + } + ], + "source": [ + "for idx in range(y_vec.shape[0]):\n", + " i = ceil(np.sqrt(2 * (idx + 1) + 0.25) - 0.5)\n", + " j = int(idx - 0.5 * i * (i - 1)) \n", + " print(y_vec[idx], y[i, j], (i, j))\n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": 76, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[ 0, 1, 2, 4, 7],\n", + " [ 1, 0, 3, 5, 8],\n", + " [ 2, 3, 0, 6, 9],\n", + " [ 4, 5, 6, 0, 10],\n", + " [ 7, 8, 9, 10, 0]])" + ] + }, + "execution_count": 76, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "y" + ] + }, + { + "cell_type": "code", + "execution_count": 70, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "1" + ] + }, + "execution_count": 70, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "y[i, j]" + ] + }, + { + "cell_type": "code", + "execution_count": 71, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "1" + ] + }, + "execution_count": 71, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "y_vec[idx]" + ] + }, + { + "cell_type": "code", + "execution_count": 35, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "1.0" + ] + }, + "execution_count": 35, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "i * (i - 1) / 2 + j" + ] + }, + { + "cell_type": "code", + "execution_count": 38, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[ 0, 1, 2, 3, 4],\n", + " [ 5, 6, 7, 8, 9],\n", + " [10, 11, 12, 13, 14],\n", + " [15, 16, 17, 18, 19],\n", + " [20, 21, 22, 23, 24]])" + ] + }, + "execution_count": 38, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "y" + ] + }, + { + "cell_type": "code", + "execution_count": 77, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([26, 2, 2, 7, 37, 41, 20, 53, 47, 1, -1, -1, -1, -1, -1, -1, -1,\n", + " -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,\n", + " -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,\n", + " -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1])" + ] + }, + "execution_count": 77, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "pos_subsample[0, 0]" + ] + }, + { + "cell_type": "code", + "execution_count": 80, + "metadata": {}, + "outputs": [], + "source": [ + "dyads = set()" + ] + }, + { + "cell_type": "code", + "execution_count": 97, + "metadata": {}, + "outputs": [ + { + "ename": "TypeError", + "evalue": "unhashable type: 'list'", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mTypeError\u001b[0m Traceback (most recent call last)", + "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mdyads\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0madd\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msorted\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2\u001b[0m \u001b[0mdyads\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0madd\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msorted\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3\u001b[0m \u001b[0mdyads\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0madd\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msorted\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m2\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;31mTypeError\u001b[0m: unhashable type: 'list'" + ] + } + ], + "source": [ + "dyads.add(tuple(sorted([0, 1]))\n", + "dyads.add(sorted([2, 1]))\n", + "dyads.add(sorted([1, 2]))" + ] + }, + { + "cell_type": "code", + "execution_count": 88, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[0, 1],\n", + " [1, 2],\n", + " [2, 1]])" + ] + }, + "execution_count": 88, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "np.asarray(list(dyads))\n" + ] + }, + { + "cell_type": "code", + "execution_count": 92, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(0, 1)" + ] + }, + "execution_count": 92, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "tuple(sorted((1, 0)))" + ] + }, + { + "cell_type": "code", + "execution_count": 95, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[0, 1]" + ] + }, + "execution_count": 95, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sorted([1, 0])" + ] + }, + { + "cell_type": "code", + "execution_count": 102, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[1, 0],\n", + " [2, 1]])" + ] + }, + "execution_count": 102, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "unique_dyads = set(tuple(sorted(l, reverse=True)) for l in dyads)\n", + "np.asarray(list(unique_dyads))" + ] + }, + { + "cell_type": "code", + "execution_count": 104, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "1151712000" + ] + }, + "execution_count": 104, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "4000 * (4000 - 1) * 3 * 24" + ] + }, + { + "cell_type": "code", + "execution_count": 100, + "metadata": {}, + "outputs": [], + "source": [ + "?sorted" + ] + }, + { + "cell_type": "code", + "execution_count": 33, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "k = 0\n", + "t = 0\n", + "y_vec = Y[k, t][np.tril_indices_from(Y[k, t], k=-1)]" + ] + }, + { + "cell_type": "code", + "execution_count": 35, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(2080,)" + ] + }, + "execution_count": 35, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "y_vec.shape" + ] + }, + { + "cell_type": "code", + "execution_count": 36, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "2080.0" + ] + }, + "execution_count": 36, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "0.5 * 65 * (65 - 1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { diff --git a/notebooks/Synthetic Data (Different Layer Types).ipynb b/notebooks/Synthetic Data (Different Layer Types).ipynb index b032179..c99d0a8 100644 --- a/notebooks/Synthetic Data (Different Layer Types).ipynb +++ b/notebooks/Synthetic Data (Different Layer Types).ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 1, + "execution_count": 26, "metadata": {}, "outputs": [], "source": [ @@ -20,23 +20,11 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 27, "metadata": { "scrolled": true }, - "outputs": [ - { - "ename": "ValueError", - "evalue": "too many values to unpack (expected 4)", - "output_type": "error", - "traceback": [ - "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)", - "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 12\u001b[0m \u001b[0mlmbda\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mlmbda\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mintercept\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1.0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 13\u001b[0m \u001b[0massortative_reference\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mFalse\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 14\u001b[0;31m random_state=423)\n\u001b[0m", - "\u001b[0;31mValueError\u001b[0m: too many values to unpack (expected 4)" - ] - } - ], + "outputs": [], "source": [ "from multidynet.datasets import simple_dynamic_multilayer_network\n", "from multidynet.plots import plot_network\n", @@ -46,11 +34,12 @@ " [-2.0, 2.0],\n", " [-1., -1.]])\n", "\n", - "Y, X, lmbda, intercept = simple_dynamic_multilayer_network(\n", - " n_nodes=100, n_time_steps=10,\n", + "Y, X, lmbda, intercept, delta = simple_dynamic_multilayer_network(\n", + " n_nodes=1000, n_time_steps=10,\n", " tau_sq=4.0, sigma_sq=0.1,\n", - " lmbda=lmbda, intercept=-1.0,\n", + " lmbda=lmbda, intercept=0.0,\n", " assortative_reference=False,\n", + " include_node_effects=True,\n", " random_state=423)" ] }, @@ -3347,18 +3336,15 @@ }, { "cell_type": "code", - "execution_count": 58, + "execution_count": null, "metadata": {}, "outputs": [ { - "data": { - "text/plain": [ - "" - ] - }, - "execution_count": 58, - "metadata": {}, - "output_type": "execute_result" + "name": "stderr", + "output_type": "stream", + "text": [ + " 3%|â–Ž | 14/500 [01:51<1:07:00, 8.27s/it]" + ] } ], "source": [ @@ -3373,15 +3359,17 @@ " lambda_var_prior=10., \n", " intercept_var_prior=10.,\n", " lambda_odds_prior=2,\n", + " include_node_effects=True,\n", " tol=1e-1,\n", - " n_init=8,\n", + " n_init=1,\n", + " n_jobs=1,\n", " random_state=123)\n", "model.fit(Y_train)" ] }, { "cell_type": "code", - "execution_count": 60, + "execution_count": 14, "metadata": {}, "outputs": [ { @@ -4164,7 +4152,7 @@ { "data": { "text/html": [ - "" + "" ], "text/plain": [ "" @@ -4179,7 +4167,7 @@ "Text(0.5, 0, 'iteration #')" ] }, - "execution_count": 60, + "execution_count": 14, "metadata": {}, "output_type": "execute_result" } @@ -4191,6 +4179,89 @@ "ax.set_xlabel('iteration #')" ] }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([ 2.58899500e-01, -1.23571196e-01, 7.06914657e-01, -1.70152924e+00,\n", + " -7.84541897e-01, 1.44719618e+00, 1.24128060e+00, -5.63144731e-01,\n", + " 4.00922219e-04, -8.27674193e-01, -1.01958696e+00, 3.25779634e-01,\n", + " -1.86201604e-01, 1.60432943e-01, -3.22497012e-01, 2.19120888e+00,\n", + " 1.58739213e+00, 2.95883490e-01, 9.06702366e-01, -6.31579972e-01,\n", + " 6.77975569e-01, 1.18971568e+00, -5.05411660e-01, 7.47895518e-01,\n", + " 3.35551375e+00, 2.71052938e-01, -3.46394145e-02, -4.84400631e-01,\n", + " 1.72188179e-01, -1.58416273e+00, -7.05639449e-02, -1.28860841e-01,\n", + " -1.36987005e+00, -1.64267978e+00, 4.11786556e-01, -6.09175293e-01,\n", + " -5.75304029e-01, -6.24189362e-01, 4.85688595e-01, -1.08374608e+00,\n", + " 1.17070449e+00, 1.31227592e+00, -1.52702204e-01, -7.34518556e-01,\n", + " 1.21496675e+00, 1.34564847e-01, -1.29692976e+00, -1.27542390e+00,\n", + " -1.08244811e+00, -2.67752479e-01, 6.31368645e-01, 9.87393798e-01,\n", + " 1.04796347e-01, 8.33305423e-01, 8.35048149e-01, -2.18969394e-01,\n", + " 1.29185590e+00, -6.01445620e-01, 2.64246204e-01, -2.54798162e-02,\n", + " -2.01010782e-01, -3.26791308e-01, -1.21510697e+00, 5.14754339e-01,\n", + " 1.87880686e+00, 8.84948885e-01, 1.44228067e-01, 2.14217622e+00,\n", + " 8.56071768e-01, -7.91608523e-01, 6.70492979e-01, -4.11619054e-01,\n", + " 1.16489224e+00, -6.61446061e-01, 1.39802375e+00, 1.13772245e+00,\n", + " -6.13446599e-01, -9.54852102e-01, -6.91488426e-02, 2.39874291e+00,\n", + " -5.38204184e-01, -2.75682041e-01, -6.17622706e-01, -6.54076767e-01,\n", + " 8.01647157e-01, 4.98012097e-01, 5.41541558e-01, 3.87243309e-01,\n", + " 3.11963350e-01, -1.26038206e+00, -1.34981885e+00, 2.22988989e-01,\n", + " -1.86095381e+00, 2.19289891e-01, -2.68171365e+00, 6.37212056e-01,\n", + " 1.06294565e+00, -3.07299616e-01, 5.55787145e-01, -1.10246282e+00])" + ] + }, + "execution_count": 15, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "model.delta_[0]" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([ 0.00803449, -0.18225511, 0.56074225, -1.76970464, -0.80514591,\n", + " 1.739377 , 1.24893226, -0.53882758, 0.14871241, -0.83306017,\n", + " -0.88482437, 0.29957807, -0.17680044, 0.20060597, -0.48300683,\n", + " 2.29152166, 1.40963996, 0.30886943, 1.10397115, -0.8505312 ,\n", + " 0.70167502, 1.19278295, -0.41097917, 0.71915733, 3.45910332,\n", + " 0.40701259, -0.33850995, -0.59780589, 0.07416049, -1.62664257,\n", + " 0.0583756 , 0.01308301, -1.41614265, -1.35984878, 0.53209186,\n", + " -0.63794355, -0.60640619, -0.60581548, 0.41398322, -1.24047692,\n", + " 1.14036145, 1.40579192, -0.18253967, -0.84248392, 1.15762719,\n", + " 0.32932595, -1.35935668, -1.29840797, -0.95694818, -0.44415006,\n", + " 0.38176558, 1.07240654, 0.15428981, 0.90512104, 0.60830456,\n", + " -0.17615639, 0.91098858, -0.57039192, 0.45922055, -0.07207278,\n", + " -0.18569592, -0.20774307, -1.06973177, 0.52350358, 2.15306017,\n", + " 0.70206151, 0.18108305, 1.90585502, 0.94840086, -0.8016874 ,\n", + " 0.84653392, -0.49577029, 0.84422595, -0.74674457, 1.13155034,\n", + " 1.01209522, -0.44568097, -0.93669264, -0.02552996, 2.46266287,\n", + " -0.4786524 , -0.28142721, -0.43523487, -0.69622877, 0.78425369,\n", + " 0.45226119, 0.59642097, 0.15365606, 0.44861798, -1.12022313,\n", + " -1.36898467, 0.15443402, -1.8713966 , 0.30531492, -2.25358515,\n", + " 0.72709477, 1.00637826, -0.41331049, 0.62907374, -1.19883164])" + ] + }, + "execution_count": 17, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "delta[0]" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -4207,19 +4278,19 @@ }, { "cell_type": "code", - "execution_count": 59, + "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 1. , 1. ],\n", - " [ 0.51249774, -0.50820869],\n", - " [-1.83194491, 1.82984983],\n", - " [-1.00120164, -1.00713605]])" + " [-0.50594665, 0.51868955],\n", + " [ 1.91832648, -1.95076181],\n", + " [-1.02121132, -1.0358395 ]])" ] }, - "execution_count": 59, + "execution_count": 18, "metadata": {}, "output_type": "execute_result" } @@ -4237,16 +4308,16 @@ }, { "cell_type": "code", - "execution_count": 69, + "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "array([-1.01446848, -1.00378947, -0.93461593, -1.0283735 ])" + "array([0., 0., 0., 0.])" ] }, - "execution_count": 69, + "execution_count": 19, "metadata": {}, "output_type": "execute_result" } @@ -4264,7 +4335,7 @@ }, { "cell_type": "code", - "execution_count": 70, + "execution_count": 20, "metadata": { "scrolled": true }, @@ -4272,10 +4343,10 @@ { "data": { "text/plain": [ - "4.586119677223207" + "4.517209610091886" ] }, - "execution_count": 70, + "execution_count": 20, "metadata": {}, "output_type": "execute_result" } @@ -4293,16 +4364,16 @@ }, { "cell_type": "code", - "execution_count": 71, + "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "0.060851727919222554" + "0.06417896024085065" ] }, - "execution_count": 71, + "execution_count": 21, "metadata": {}, "output_type": "execute_result" } @@ -4320,7 +4391,7 @@ }, { "cell_type": "code", - "execution_count": 96, + "execution_count": 25, "metadata": {}, "outputs": [ { @@ -5103,7 +5174,7 @@ { "data": { "text/html": [ - "" + "" ], "text/plain": [ "" @@ -5137,7 +5208,7 @@ }, { "cell_type": "code", - "execution_count": 61, + "execution_count": 23, "metadata": { "scrolled": true }, @@ -5145,10 +5216,10 @@ { "data": { "text/plain": [ - "0.9646799483170649" + "0.9671974701938675" ] }, - "execution_count": 61, + "execution_count": 23, "metadata": {}, "output_type": "execute_result" } @@ -5166,16 +5237,16 @@ }, { "cell_type": "code", - "execution_count": 62, + "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "0.962888086908839" + "0.9654683802943135" ] }, - "execution_count": 62, + "execution_count": 22, "metadata": {}, "output_type": "execute_result" } @@ -5195,7 +5266,7 @@ }, { "cell_type": "code", - "execution_count": 63, + "execution_count": 24, "metadata": {}, "outputs": [ { @@ -5978,7 +6049,7 @@ { "data": { "text/html": [ - "" + "" ], "text/plain": [ "" @@ -5993,7 +6064,7 @@ "Text(0.5, 0, 'Estimated Probability ($\\\\hat{\\\\pi}_{ijt}^k$)')" ] }, - "execution_count": 63, + "execution_count": 24, "metadata": {}, "output_type": "execute_result" }