1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
|
# cython: boundscheck=False, cdivision=True
import numpy as np
cimport numpy as np
from libc.math cimport log, exp
DTYPE = np.float64
ctypedef np.float_t DTYPE_t
cdef DTYPE_t plogis(DTYPE_t weight, DTYPE_t delta):
return 1./(1. + exp(-weight/delta))
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
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):
"""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
result = log(1. - structural + structural * temporal)
# print 'stnv', structural, temporal
return result
def ml3(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
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 = 148152
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] probs_nv = np.zeros(len(non_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)
# loop through victims
for i, parents in enumerate(victims.itervalues()):
# 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)
for (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]
dists = [dist for (dist, dt, w1, w2, w3) in parents]
dts = [dt for (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]
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]
# 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_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.09, 1., .05):
thresh = log(beta/(3012.*(1.-beta)))
# print 'beta:', beta, 'thresh:', thresh, 'infected:', len(probs[probs>=thresh])
roots = n_roots + len(probs[probs<thresh])
beta_add = 0.
# add probability for realized edges and subtract probability these edges fail
beta_add += (probs[probs>=thresh]).sum()
# add probability for the seeds and non-seeds
beta_add += roots * log(beta) + len(probs[probs>=thresh]) * log(1. - beta)
if beta_add > max_beta_add:
max_beta = beta
max_roots = roots
max_beta_add = beta_add
# 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
return (beta, roots, ll)
|
