Coverage for Python_files/predictor_tools.py: 87%

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

39 statements  

1""" 

2This code is aimed to provide tools for prediction process.  

3""" 

4 

5import numpy as np 

6 

7 

8 

9def predictions(params, history, alpha, mu, T = None): 

10 """ 

11 Returns the expected total numbers of points for a set of time points 

12  

13 params -- parameter tuple (p,beta) of the Hawkes process 

14 history -- (n,2) numpy array containing marked time points (t_i,m_i)  

15 alpha -- power parameter of the power-law mark distribution 

16 mu -- min value parameter of the power-law mark distribution 

17 T -- 1D-array of times (i.e ends of observation window) 

18 """ 

19 if T : 

20 if not isinstance(T, (float, int) ) or T < 0: 

21 raise Exception(" T must be an float or int greater than 0") 

22 

23 if not isinstance(params, np.ndarray) : 

24 raise Exception(" params must be a np.ndarray") 

25 

26 if not isinstance(params[0], np.floating) : 

27 raise Exception ("p must be a int or float") 

28 

29 if not isinstance(params[0], np.floating) or params[1]<0: 

30 raise Exception ("beta must be a int or float greater than 0") 

31 

32 if not isinstance(history, np.ndarray) or history.shape[1]!=2 : 

33 raise Exception(" history must be an np.array with following shape : (n,2)") 

34 

35 if not isinstance(alpha, (int,float)): 

36 raise Exception(" alpha must be an float or int ") 

37 

38 if not isinstance(mu, (int,float)): 

39 raise Exception(" mu must be an float or int ") 

40 p,beta = params 

41 

42 tis = history[:,0] 

43 if T is None: 

44 T = np.linspace(60,tis[-1],1000) 

45 

46 N = np.zeros((len(T),2)) 

47 N[:,0] = T 

48 

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

50 n_star = p * EM 

51 if n_star >= 1: 

52 raise Exception(f"Branching factor {n_star:.2f} greater than one") 

53 

54 Si, ti_prev, i = 0., 0., 0 

55 

56 for j,t in enumerate(T): 

57 for (ti,mi) in history[i:]: 

58 if ti >= t: 

59 break 

60 else: 

61 Si = Si * np.exp(-beta * (ti - ti_prev)) + mi 

62 ti_prev = ti 

63 i += 1 

64 

65 n = i + 1 

66 G1 = p * Si * np.exp(-beta * (t - ti_prev)) 

67 N[j,1] = n + G1 / (1. - n_star) 

68 return N,n_star,G1