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

84 statements  

1""" this is the class version of the previous code""" 

2 

3import numpy as np 

4import scipy.optimize as optim 

5 

6class Cascade : 

7 def __init__(self,json) : 

8 

9 ############################################## 

10 ########### Given Attributes ############## 

11 ############################################## 

12 

13 self._json=json.value 

14 self._key = json.value["key"] 

15 self._history=json.value["tweets"] 

16 self._t=json.value["T_obs"] 

17 

18 ############################################## 

19 ##### Constants given by Mishra et al ###### 

20 ############################################## 

21 

22 

23 self._mu=1 

24 self._alpha=2.016 

25 

26 ############################################## 

27 ############ Default constants ############ 

28 ############################################## 

29 

30 

31 self._p = 0 # default value 

32 self._beta= 0 #default value 

33 self._N=1 #default value  

34 

35 

36 

37 ################################# 

38 ########## Accessors ######### 

39 ################################# 

40 

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 

55 

56 

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 

70 

71 ################################# 

72 ########## Methods ########### 

73 ################################# 

74 

75 

76 def loglikelihood(self): 

77 """ 

78 Returns the loglikelihood of a Hawkes process with exponential kernel 

79 computed with a linear time complexity 

80  

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 """ 

85 

86 p,beta,history,t = self.Get_p(),self.Get_beta(),self.Get_history(),self.Get_T() 

87 

88 

89 if p <= 0 or p >= 1 or beta <= 0.: return -np.inf 

90 

91 n = len(history) 

92 tis = history[:,0] 

93 mis = history[:,1] 

94 

95 LL = (n-1) * np.log(p * beta) 

96 logA = -np.inf 

97 prev_ti, prev_mi = history[0] 

98 

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)) 

103 

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 

108 

109 logA = np.log(prev_mi + np.exp(logA)) - beta * (t - prev_ti) 

110 LL -= p * (np.sum(mis) - np.exp(logA)) 

111 

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) 

118 

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]]) 

140 

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. 

144 

145 # Compute the covariance inverse (precision matrix) once for all 

146 inv_cov_prior = np.asmatrix(cov_prior).I 

147 

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) 

151 

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) 

161 

162 EM = mu * (alpha - 1) / (alpha - 2) 

163 eps = 1.E-8 

164 

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) 

168 

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 ) 

174 

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