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

36 statements  

« prev     ^ index     » next       coverage.py v7.11.0, created at 2026-05-15 16:30 +0000

1""" 

2""" 

3import numpy as np 

4from causalspyne.implicit_omega_condition_on_w import gen_joint_w_omega, get_max_degree 

5from causalspyne.implicit_omega import get_extreme_eigenvalue 

6 

7 

8def cov_to_corr(covariance_matrix): 

9 # Calculate the standard deviations 

10 std_devs = np.sqrt(np.diag(covariance_matrix)) 

11 

12 # Create a diagonal matrix of inverse standard deviations 

13 inv_std_devs = np.diag(1 / std_devs) 

14 

15 # Calculate the correlation matrix 

16 correlation_matrix = inv_std_devs @ covariance_matrix @ inv_std_devs 

17 

18 return correlation_matrix 

19 

20 

21def congruent(mat_op, mat): 

22 """ 

23 congruent transform is defined to be (mat_op)^T mat (mat_op) 

24 """ 

25 first2 = np.matmul(mat_op.T, mat) 

26 mat_rst = np.matmul(first2, mat_op) 

27 return mat_rst 

28 

29 

30def gen_sigma_y(max_w=0.5): 

31 """ 

32 generate Sigma matrix according to $${(I-W)}^{-1}\Omega{(I-W)}^{-1}$$ 

33 """ 

34 mat_w_binary, mat_omega = gen_joint_w_omega() 

35 mat_weight = np.random.uniform( 

36 low=-max_w, high=max_w, size=mat_w_binary.shape) 

37 

38 mat_w = mat_w_binary * mat_weight 

39 print(f"adjacency weighted: {mat_w}") 

40 

41 dmax = get_max_degree(mat_w) 

42 kernel = np.identity(mat_w.shape[0]) - mat_w 

43 inv_id_minus_w = np.linalg.inv(kernel) 

44 mat_sigma = congruent(inv_id_minus_w.T, mat_omega) 

45 max_omega = np.max(mat_omega) 

46 return mat_sigma, dmax, max_omega 

47 

48 

49def gen_spectrum(): 

50 mat_sigma, dmax, max_omega = gen_sigma_y() 

51 print(f"sigma: {mat_sigma}") 

52 inv_sigma = np.linalg.inv(mat_sigma) 

53 eigv_max, _ = get_extreme_eigenvalue(inv_sigma) 

54 print(f"max eigen value: {eigv_max}") 

55 print(f"max degree: {dmax}, max_omega: {max_omega}") 

56 mat_precision = cov_to_corr(inv_sigma) 

57 eigv_max, _ = get_extreme_eigenvalue(mat_precision) 

58 print(f"max eigen value: {eigv_max}") 

59 return eigv_max 

60 

61 

62if __name__ == "__main__": 

63 gen_spectrum()