in which the model works!!

1 files changed, 27 insertions, 29 deletions

diff --git a/experiments/ml.pyx b/experiments/ml.pyx
index c8c99cc..858207f 100644
--- a/experiments/ml.pyx
+++ b/experiments/ml.pyx
@@ -10,14 +10,9 @@ cdef DTYPE_t weight_success(int dist, int dt, DTYPE_t alpha, DTYPE_t delta,
                             DTYPE_t w1, DTYPE_t w2, DTYPE_t w3):
     """weight for successful infection, exponential time model"""
     cdef DTYPE_t structural, temporal, result
-    structural = dist * log(delta)
-    # structural = plogis(w1,delta) * plogis(w2,delta) * plogis(w3,delta)
-    temporal = log(exp(alpha)-1.) - alpha*dt
-    # temporal = 1 - exp(-alpha*dt)
-    # if exp(-alpha*dt)==0.: print 'UNDERFLOW ERROR'
-    # temporal = 1. / (1. + (dt - 1.)/alpha)**0.01 - 1. / (1. + dt/alpha)**0.01
-    result = structural + temporal
-    # print 'st', structural, temporal
+    structural = delta ** dist
+    temporal = log(exp(alpha)-1.) - alpha*dt/7.
+    result = log(structural) + temporal
     return result
 
 cdef DTYPE_t weight_failure(int dist, int dt, DTYPE_t alpha, DTYPE_t delta,
@@ -25,17 +20,14 @@ cdef DTYPE_t weight_failure(int dist, int dt, DTYPE_t alpha, DTYPE_t delta,
     """weight for failed infection, exponential time model"""
     cdef DTYPE_t structural, temporal, result
     structural = delta ** dist
-    # structural = plogis(w1,delta) * plogis(w2,delta) * plogis(w3,delta)
-    temporal = exp(-alpha * dt)
-    # temporal = 1.  - 1. / (1. + dt/alpha)**0.01
+    temporal = exp(-alpha * dt/7.)
     result = log(1. - structural + structural * temporal)
-    # print 'stnv', structural, temporal
     return result
 
 def ml(dict root_victims, dict victims, dict non_victims, DTYPE_t age,
        DTYPE_t alpha, DTYPE_t delta):
     cdef:
-        int n_roots, n_victims, n_nodes, roots, i, dist, dt, t, l
+        int n_roots, n_victims, roots, i, dist, dt, t, l
         DTYPE_t beta, ll, beta_add, max_beta, max_beta_add
         list parents, failures, successes
     n_roots, n_victims = len(root_victims), len(victims)
@@ -43,6 +35,8 @@ def ml(dict root_victims, dict victims, dict non_victims, DTYPE_t age,
     cdef:
         np.ndarray[DTYPE_t] probs = np.zeros(n_victims, dtype=DTYPE)
         np.ndarray[DTYPE_t] probs_fail = np.zeros(n_victims, dtype=DTYPE)
+        np.ndarray[DTYPE_t] parent_dists = np.zeros(n_victims, dtype=DTYPE)
+        np.ndarray[DTYPE_t] parent_dts = np.zeros(n_victims, dtype=DTYPE)
         np.ndarray[DTYPE_t] probs_nv = np.zeros(len(non_victims), dtype=DTYPE)
 
     # loop through victims
@@ -51,43 +45,47 @@ def ml(dict root_victims, dict victims, dict non_victims, DTYPE_t age,
         # 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)
-                    for (dist, dt, w1, w2, w3) in parents]
+                    for (prnt, dist, dt, w1, w2, w3) in parents]
         probs_fail[i] = sum(failures)
         successes = [weight_success(dist, dt, alpha, delta, w1, w2, w3)
-                     for (dist, dt, w1, w2, w3) in parents]
+                     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] = max(s - failures[l] for l, s in enumerate(successes)) 
+        probs[i] = float("-inf")
+        for l, s in enumerate(successes):
+            prob = s - failures[l]
+            if prob > probs[i]:
+                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)
-                    for (dist, dt, w1, w2, w3) in parents]
+                    for (prnt, dist, dt, w1, w2, w3) in parents]
         probs_nv[i] = sum(failures)
 
-    # print successes
-    # print failures
-    # print probs
-
     # calculate log likelihood
-    # probs.sort(); probs = probs[::-1] # sort probs in descending order
-    # cdef:
-    #     np.ndarray[DTYPE_t] cums = probs.cumsum()
-    ll =  probs_fail.sum() # add probability that all edges to victims fail
+    ll =  probs_fail.sum() # add probability that all edges to all victims fail
     ll += probs_nv.sum() # add probability that all edges to non_victims fail
 
     # print 'probs', probs
     max_beta_add = float('-inf')
     # iterate over all victim nodes to find the optimal threshold
-    for beta in np.arange(0.001, .2, .002):
+    for beta in np.arange(0.01, 1., .01):
         thresh = log(beta/(3012.*(1.-beta)))
-        # print 'beta:', beta, 'thresh:', thresh, 'infected:', len(probs[probs>=thresh])
-        roots = n_roots + len(probs[probs<thresh])
+        seeds = probs<thresh
+        non_seeds = probs>=thresh
+        roots = n_roots + sum(seeds)
 
         beta_add = 0.
         # add probability for realized edges and subtract probability these edges fail
-        beta_add += (probs[probs>=thresh]).sum()
+        beta_add += (probs[non_seeds]).sum()
         # add probability for the seeds and non-seeds
         beta_add += roots * log(beta/3012.) + (n_nodes-roots) * log(1. - beta)