Coverage for Python_files/cascade_class.py: 0%
Shortcuts on this page
r m x toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
Shortcuts on this page
r m x toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1""" this is the class version of the previous code"""
3import numpy as np
4import scipy.optimize as optim
6class Cascade :
7 def __init__(self,json) :
9 ##############################################
10 ########### Given Attributes ##############
11 ##############################################
13 self._json=json.value
14 self._key = json.value["key"]
15 self._history=json.value["tweets"]
16 self._t=json.value["T_obs"]
18 ##############################################
19 ##### Constants given by Mishra et al ######
20 ##############################################
23 self._mu=1
24 self._alpha=2.016
26 ##############################################
27 ############ Default constants ############
28 ##############################################
31 self._p = 0 # default value
32 self._beta= 0 #default value
33 self._N=1 #default value
37 #################################
38 ########## Accessors #########
39 #################################
41 def Get_p(self):
42 return self._p
43 def Get_beta(self):
44 return self._beta
45 def Get_N(self):
46 return self._N
47 def Get_history(self):
48 return self._history
49 def Get_T(self):
50 return self._t
51 def Get_alpha(self):
52 return self._alpha
53 def Get_mu(self):
54 return self._mu
57 #################################
58 ########## Setters ########
59 #################################
60 def Set_p(self,value):
61 self._p=value
62 def Set_beta(self,value):
63 self._beta=value
64 def Set_N(self,value):
65 self._N=value
66 def Set_history(self,value):
67 self._history=value
68 def Set_T(self,value):
69 self._t=value
71 #################################
72 ########## Methods ###########
73 #################################
76 def loglikelihood(self):
77 """
78 Returns the loglikelihood of a Hawkes process with exponential kernel
79 computed with a linear time complexity
81 params -- parameter tuple (p,beta) of the Hawkes process
82 history -- (n,2) numpy array containing marked time points (t_i,m_i)
83 t -- current time (i.e end of observation window)
84 """
86 p,beta,history,t = self.Get_p(),self.Get_beta(),self.Get_history(),self.Get_T()
89 if p <= 0 or p >= 1 or beta <= 0.: return -np.inf
91 n = len(history)
92 tis = history[:,0]
93 mis = history[:,1]
95 LL = (n-1) * np.log(p * beta)
96 logA = -np.inf
97 prev_ti, prev_mi = history[0]
99 i = 0
100 for ti,mi in history[1:]:
101 if(prev_mi + np.exp(logA) <= 0):
102 print("Bad value", prev_mi + np.exp(logA))
104 logA = np.log(prev_mi + np.exp(logA)) - beta * (ti - prev_ti)
105 LL += logA
106 prev_ti,prev_mi = ti,mi
107 i += 1
109 logA = np.log(prev_mi + np.exp(logA)) - beta * (t - prev_ti)
110 LL -= p * (np.sum(mis) - np.exp(logA))
112 return LL
113 def compute_MAP(self,
114 prior_params = [ 0.02, 0.0002, 0.01, 0.001, -0.1],
115 max_n_star = 1, display=False):
116 """
117 Returns the pair of the estimated logdensity of a posteriori and parameters (as a numpy array)
119 history -- (n,2) numpy array containing marked time points (t_i,m_i)
120 t -- current time (i.e end of observation window)
121 alpha -- power parameter of the power-law mark distribution
122 mu -- min value parameter of the power-law mark distribution
123 prior_params -- list (mu_p, mu_beta, sig_p, sig_beta, corr) of hyper parameters of the prior
124 -- where:
125 -- mu_p: is the prior mean value of p
126 -- mu_beta: is the prior mean value of beta
127 -- sig_p: is the prior standard deviation of p
128 -- sig_beta: is the prior standard deviation of beta
129 -- corr: is the correlation coefficient between p and beta
130 max_n_star -- maximum authorized value of the branching factor (defines the upper bound of p)
131 display -- verbose flag to display optimization iterations (see 'disp' options of optim.optimize)
132 """
133 history,t = self.Get_history(),self.Get_T()
134 alpha,mu=self.Get_alpha(),self.Get_mu()
135 # Compute prior moments
136 mu_p, mu_beta, sig_p, sig_beta, corr = prior_params
137 sample_mean = np.array([mu_p, mu_beta])
138 cov_p_beta = corr * sig_p * sig_beta
139 Q = np.array([[sig_p ** 2, cov_p_beta], [cov_p_beta, sig_beta **2]])
141 # Apply method of moments
142 cov_prior = np.log(Q / sample_mean.reshape((-1,1)) / sample_mean.reshape((1,-1)) + 1)
143 mean_prior = np.log(sample_mean) - np.diag(cov_prior) / 2.
145 # Compute the covariance inverse (precision matrix) once for all
146 inv_cov_prior = np.asmatrix(cov_prior).I
148 # Define the target function to minimize as minus the log of the a posteriori density
149 def target(params):
150 log_params = np.log(params)
152 if np.any(np.isnan(log_params)):
153 return np.inf
154 else:
155 dparams = np.asmatrix(log_params - mean_prior)
156 prior_term = float(- 1/2 * dparams * inv_cov_prior * dparams.T)
157 self.Set_beta(params[1])
158 self.Set_p(params[0])
159 logLL = self.loglikelihood()
160 return - (prior_term + logLL)
162 EM = mu * (alpha - 1) / (alpha - 2)
163 eps = 1.E-8
165 # Set realistic bounds on p and beta
166 p_min, p_max = eps, max_n_star/EM - eps
167 beta_min, beta_max = 1/(3600. * 24 * 10), 1/(60. * 1)
169 # Define the bounds on p (first column) and beta (second column)
170 bounds = optim.Bounds(
171 np.array([p_min, beta_min]),
172 np.array([p_max, beta_max])
173 )
175 # Run the optimization
176 res = optim.minimize(
177 target, sample_mean,
178 method='Powell',
179 bounds=bounds,
180 options={'xtol': 1e-8, 'disp': display}
181 )
182 # Returns the loglikelihood and found parameters
183 return(-res.fun, res.x) ## res.x contains p et beta