Coverage for src/causalspyne/implicit_omega.py: 0%

1""" 

2This module is used to generate the implicit omega matrix. 

3 

4To generate a symmetric matrix, 

5 - A+A^T 

6 - np.tril(A) + np.tril(A, -1).T 

7 np.tril(A, -1) returns a copy of the input array A with elements above the 

8 first sub-diagonal zeroed out 

9""" 

10import numpy as np 

11 

12 

13def is_diagonal(matrix): 

14 return np.count_nonzero(matrix - np.diag(np.diag(matrix))) == 0 

15 

16 

17def is_diag_dom(mat): 

18 for i_row, row in enumerate(mat): 

19 energy_diag = np.abs(mat[i_row, i_row]) 

20 energy_row_off_diag = np.sum(np.abs(row)) - energy_diag 

21 flag = energy_diag <= energy_row_off_diag 

22 if flag: 

23 return False 

24 return True 

25 

26 

27def gen_omega(n_ob_v, max_omega, delta=1.0, p_sparse=0.1): 

28 """ 

29 :n_ob_v: number of observed variables 

30 :max_omega: maximum value of omega 

31 :returns: None 

32 """ 

33 assert max_omega > 0 

34 # Generate Bernoulli random numbers 

35 mat_bernoulli = np.random.binomial( 

36 n=1, p=p_sparse, size=(n_ob_v, n_ob_v)) 

37 mat_omega = np.random.uniform( 

38 low=0.1, high=max_omega, size=(n_ob_v, n_ob_v)) 

39 mat_sparse_omega = mat_bernoulli * mat_omega 

40 mat_sym = 0.5 * (mat_sparse_omega + mat_sparse_omega.T) 

41 

42 if is_diagonal(mat_sym): 

43 rng = np.random.default_rng() 

44 sampled_row_ind = rng.integers(low=0, high=n_ob_v - 1, size=1) 

45 off_diag_v = mat_omega[sampled_row_ind, sampled_row_ind + 1] 

46 mat_sym[sampled_row_ind, sampled_row_ind + 1] = off_diag_v 

47 mat_sym[sampled_row_ind + 1, sampled_row_ind] = off_diag_v 

48 for i_row, row in enumerate(mat_sym): 

49 old_diag_v = np.abs(mat_sym[i_row, i_row]) 

50 energy_row_off_diag = np.sum(np.abs(row)) - np.abs(old_diag_v) 

51 mat_sym[i_row, i_row] = delta + energy_row_off_diag 

52 assert is_diag_dom(mat_sym) 

53 return mat_sym, delta + max_omega 

54 

55 

56def get_extreme_eigenvalue(matrix): 

57 """ 

58 return: max, min 

59 """ 

60 eigenvalues = np.linalg.eigvals(matrix) 

61 return np.max(np.real(eigenvalues)), np.min(np.real(eigenvalues)) 

62 

63 

64if __name__ == "__main__": 

65 max_omega = 0.7 

66 mat, rho = gen_omega(3, max_omega) 

67 print(mat) 

68 max_eig, min_eig = get_extreme_eigenvalue(mat) 

69 print(f"min, max eigenvalue {min_eig, max_eig}") 

70 print(rho, 2 * rho)