updated temporal model for different time scales, ran some tests

1 files changed, 26 insertions, 19 deletions

diff --git a/experiments/ml.pyx b/experiments/ml.pyx
index 0d72c44..7eca30a 100644
--- a/experiments/ml.pyx
+++ b/experiments/ml.pyx
@@ -6,32 +6,38 @@ from libc.math cimport log, exp
 DTYPE = np.float64
 ctypedef np.float_t DTYPE_t
 
-cdef DTYPE_t weight_success(int dist, int dt, DTYPE_t alpha, DTYPE_t delta,
-                            DTYPE_t w1, DTYPE_t w2, DTYPE_t w3):
+cdef DTYPE_t weight_success(int dist, int dt, DTYPE_t alpha, DTYPE_t delta, DTYPE_t lmbda,
+                            DTYPE_t t_scale, DTYPE_t w1, DTYPE_t w2, DTYPE_t w3):
     """weight for successful infection, exponential time model"""
     cdef DTYPE_t structural, temporal, result
     structural = delta ** dist
-    temporal = log(exp(alpha)-1.) - alpha*dt/1.
+    # structural = delta/(1. + 1./(w1*lmbda) + 1./(w2*lmbda) + 1./(w3*lmbda))
+    temporal = log(exp(alpha/t_scale)-1.) - alpha*dt/t_scale
+    # temporal = 1. / (1. + (dt - 1.)/alpha)**0.01 - 1. / (1. + dt/alpha)**0.01
     result = log(structural) + temporal
     return result
 
-cdef DTYPE_t weight_failure(int dist, int dt, DTYPE_t alpha, DTYPE_t delta,
-                            DTYPE_t w1, DTYPE_t w2, DTYPE_t w3):
+cdef DTYPE_t weight_failure(int dist, int dt, DTYPE_t alpha, DTYPE_t delta, DTYPE_t lmbda,
+                            DTYPE_t t_scale, DTYPE_t w1, DTYPE_t w2, DTYPE_t w3):
     """weight for failed infection, exponential time model"""
     cdef DTYPE_t structural, temporal, result
     structural = delta ** dist
-    temporal = exp(-alpha * dt/1.)
+    # structural = delta/(1. + 1./(w1*lmbda) + 1./(w2*lmbda) + 1./(w3*lmbda))
+    temporal = exp(-alpha * dt/t_scale)
+    # temporal = 1.  - 1. / (1. + dt/alpha)**0.01
     result = log(1. - structural + structural * temporal)
     return result
 
-def ml(dict root_victims, dict victims, dict non_victims, DTYPE_t age,
-       DTYPE_t alpha, DTYPE_t delta):
+def ml(dict root_victims, dict victims, dict non_victims,
+       DTYPE_t alpha, DTYPE_t delta, DTYPE_t lmbda):
     cdef:
-        int n_roots, n_victims, roots, i, dist, dt, t, l
+        int n_roots, n_victims, roots, i, dist, dt, t, l, n_days
         DTYPE_t beta, ll, beta_add, max_beta, max_beta_add
         list parents, failures, successes
     n_roots, n_victims = len(root_victims), len(victims)
-    n_nodes = 5000
+    n_nodes = 148152
+    n_days = 3012
+    t_scale = 1.
     cdef:
         np.ndarray[DTYPE_t] probs = np.zeros(n_victims, dtype=DTYPE)
         np.ndarray[DTYPE_t] probs_fail = np.zeros(n_victims, dtype=DTYPE)
@@ -44,15 +50,14 @@ def ml(dict root_victims, dict victims, dict non_victims, DTYPE_t age,
         # for each victim node i, compute the probability that all its parents
         # fail to infect it, also computes the probability that its most
         # likely parent infects it
-        failures = [weight_failure(dist, dt, alpha, delta, w1, w2, w3)
+        failures = [weight_failure(dist, dt, alpha, delta, lmbda, t_scale, w1, w2, w3)
                     for (prnt, dist, dt, w1, w2, w3) in parents]
         probs_fail[i] = sum(failures)
-        successes = [weight_success(dist, dt, alpha, delta, w1, w2, w3)
+        successes = [weight_success(dist, dt, alpha, delta, lmbda, t_scale, w1, w2, w3)
                      for (prnt, dist, dt, w1, w2, w3) in parents]
         dists = [dist for (prnt, dist, dt, w1, w2, w3) in parents]
         dts = [dt for (prnt, dist, dt, w1, w2, w3) in parents]
         # find parent that maximizes log(p) - log(\tilde{p})
-        # probs[i] = max(s - failures[l] for l, s in enumerate(successes)) 
         probs[i] = float("-inf")
         for l, s in enumerate(successes):
             prob = s - failures[l]
@@ -60,13 +65,12 @@ def ml(dict root_victims, dict victims, dict non_victims, DTYPE_t age,
                 probs[i] = prob
                 parent_dists[i] = dists[l]
                 parent_dts[i] = dts[l]
-                # probs_fail[i] = failures[l]
 
     # loop through non-victims
     for i, parents in enumerate(non_victims.itervalues()):
         # for each non victim node, compute the probability that all its
         # parents fail to infect it
-        failures = [weight_failure(dist, dt, alpha, delta, w1, w2, w3)
+        failures = [weight_failure(dist, dt, alpha, delta, lmbda, t_scale, w1, w2, w3)
                     for (prnt, dist, dt, w1, w2, w3) in parents]
         probs_nv[i] = sum(failures)
 
@@ -76,8 +80,8 @@ def ml(dict root_victims, dict victims, dict non_victims, DTYPE_t age,
 
     max_beta_add = float('-inf')
     # iterate over all victim nodes to find the optimal threshold
-    for beta in np.arange(0.01, 1., .01):
-        thresh = log(beta/(1000.*(1.-beta)))
+    for beta in np.arange(0.01, 1., 0.01):
+        thresh = log(beta/(n_days*(1.-beta)))
         seeds = probs<thresh
         non_seeds = probs>=thresh
         roots = n_roots + sum(seeds)
@@ -86,16 +90,19 @@ def ml(dict root_victims, dict victims, dict non_victims, DTYPE_t age,
         # add probability for realized edges and subtract probability these edges fail
         beta_add += (probs[non_seeds]).sum()
         # add probability for the seeds and non-seeds
-        beta_add += roots * log(beta/1000.) + (n_nodes-roots) * log(1. - beta)
+        beta_add += roots * log(beta/n_days) + (n_nodes-roots) * log(1. - beta)
 
         if beta_add > max_beta_add:
             max_beta = beta
             max_roots = roots
             max_beta_add = beta_add
+            pdists = (parent_dists[non_seeds]).mean()
+            pdts = (parent_dts[non_seeds]).mean()
             # print 'beta:', max_beta, 'add:', max_beta_add, 'roots:', max_roots
 
     ll += max_beta_add
     roots = max_roots
     beta = max_beta
-    # print n_nodes, n_roots, n_victims, max_i, roots
+    print 'dist:', pdists
+    print 'dt:', pdts
     return (beta, roots, ll)