<jupyter_start><jupyter_text>Causal Discovery under Causal SufficiencyIn certain case, it is reasonable to assume that there are no latent confounders and no selection bias.In this case, we assume that given a dataset, the true underlying causal diagram, over the measured variables,is a directed acyclic graph (DAG).A common algorithm used in these cases is the PC algorithm (named after its inventors <NAME> and <NAME>).An efficient algorithm is the recursive autonomy identification (RAI) algorithm.In the large sample limit, both algorithms a proved to recover the equivalence class of the true underlying DAG.This equivalence class is called completed partially directed graph (CPDAG), or essential graph.For demonstrating PC and RAI, we follow these steps. Initially, the required classes and methods are imported.<jupyter_code>import sys
sys.path.append('..')

import numpy as np
from causal_discovery_utils.cond_indep_tests import CondIndepCMI
from causal_discovery_algs import LearnStructRAI, LearnStructPC
from causal_discovery_utils.data_utils import get_var_size
from graphical_models import DAG, PDAG
from causal_discovery_utils.performance_measures import structural_hamming_distance_cpdag, score_bdeu
from experiment_utils.threshold_select_ci_test import search_threshold_bdeu
from matplotlib import pyplot as plt<jupyter_output><empty_output><jupyter_text>Experiment SetupFirstly, we will use the ALARM monitoring system (Beinlich et al., 1989) as the true underlying graph.It is a Bayesian network consisting of 37 nodes (8 diagnoses, 16 findings and 13 intermediate variables), and 46 edges.As an example, we use one dataset sampled from this Bayesian network by Tsamardinos et al. (2006).The full database they created, covering a range of Bayesian network and a range of datasets size can be found inthe [supplementary material for the max-min hill-climbing (MMHC) algorithm](https://pages.mtu.edu/~lebrown/supplements/mmhc_paper/mmhc_index.html).  InitializationInitially, we define the location and file name of the datasets and graph structure.<jupyter_code>data_name = 'Alarm1'  # the ALARM monitoring network (Beinlich et at., 1989)
data_size = 500  # dataset size for learning
data_ver = 1  # dataset version

data_folder = '../example_data/' + data_name + '_data/'
test_data_folder = '../example_data/' + data_name + '_testdata/'<jupyter_output><empty_output><jupyter_text>Load Training DataThe training data will be used to learn the underlying graph structure.<jupyter_code>data_train = np.loadtxt(data_folder + data_name + '_s' + str(data_size) + '_v' + str(data_ver) + '.txt',
                        dtype=int)
n_samples, n_vars = data_train.shape  # data is assumed a numpy 2d-array
graph_nodes = set(range(n_vars))  # create a set containing the nodes indices<jupyter_output><empty_output><jupyter_text>The loaded data, `data_train`, is a 2D numpy array, where its first axis is the sample indexand its second axis is the variable index.  Learn the Equivalence Class (CPDAG) of the Underlying DAGsInitially, we set the conditional independence test.We select the conditional mutual information (CMI) as it is suitable for estimating the level of correlationbetween discrete variables (all the variables in the ALARM network are discrete).<jupyter_code>CITest = CondIndepCMI  # class of the ci test to be used<jupyter_output><empty_output><jupyter_text>In addition, we search the CMI threshold that maximizes the likelihood of the graph.This strategy was suggested by Yehezkel and Lerner (2009) for the RAI algorithm. We use it herefor the PC algorithm as well. The method that searches for the threshold is `search_threshold_bdeu`.It utilizes the caching mechanism of the CI test class to search the threshold efficiently.Set the list of candidate thresholds.<jupyter_code>th_range = [i / 10000 + 0.01 for i in range(100)]  # list of candidate thresholds for CI test<jupyter_output><empty_output><jupyter_text>Find the threshold for RAI<jupyter_code>th_rai, all_scores_rai = search_threshold_bdeu(LearnStructRAI, data_train, CITest, th_range)
print('Selected RAI threshold = {:.4f}'.format(th_rai))<jupyter_output>Selected RAI threshold = 0.0159<jupyter_text>Find the threshold for PC<jupyter_code>th_pc, all_scores_pc = search_threshold_bdeu(LearnStructPC, data_train, CITest, th_range)
print('Selected PC threshold = {:.4f}'.format(th_pc))<jupyter_output>Selected PC threshold = 0.0120<jupyter_text>Learn using the RAI algorithmFirst, instantiate a CI test with the selected threshold and training dataset<jupyter_code>ci_test_rai = CITest(dataset=data_train, threshold=th_rai, count_tests=True)  # conditional independence test<jupyter_output><empty_output><jupyter_text>Then, instantiate an RAI learner<jupyter_code>rai = LearnStructRAI(nodes_set=graph_nodes, ci_test=ci_test_rai)<jupyter_output><empty_output><jupyter_text>Finally, learn the structure<jupyter_code>rai.learn_structure()  # learn structure<jupyter_output><empty_output><jupyter_text>Learn using the PC algorithmInstantiate a CI test and a PC learner, and learn the structure.<jupyter_code>ci_test_pc = CITest(dataset=data_train, threshold=th_pc, count_tests=True)  # conditional independence test
pc = LearnStructPC(nodes_set=graph_nodes, ci_test=ci_test_pc)
pc.learn_structure()  # learn structure<jupyter_output><empty_output><jupyter_text>Examine ResultsWe compare the PC and RAI algorithms using three measures:1. *Complexity* in terms of the number of CI tests required for learning the graph2. *Structural Hamming distance* of the learned graph from the true underlying CPDAG3. *BDeu* score of the learned graph calculated using a novel, large, test datasetInitially, load the true underlying DAG and calculate its equivalence class CPDAG.<jupyter_code># Load True Graph
true_dag_np = np.loadtxt(data_folder + data_name + '_graph.txt', dtype=np.int)
true_dag = DAG(graph_nodes)
true_dag.init_from_adj_mat(true_dag_np)
true_cpdag = PDAG(true_dag.nodes_set)
true_dag.convert_to_cpdag(true_cpdag)  # create equivalence class of th true DAG

# Load Test Data
data_test = np.loadtxt(
    test_data_folder + data_name + '_testdata_s' + str(data_size) + '_v' + str(data_ver) + '.txt',
    dtype=int)<jupyter_output>/var/folders/68/h1wrnzx56gvdj8f4mm8npqgr0000gp/T/ipykernel_67812/2780770500.py:2: DeprecationWarning: `np.int` is a deprecated alias for the builtin `int`. To silence this warning, use `int` by itself. Doing this will not modify any behavior and is safe. When replacing `np.int`, you may wish to use e.g. `np.int64` or `np.int32` to specify the precision. If you wish to review your current use, check the release note link for additional information.
Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations
  true_dag_np = np.loadtxt(data_folder + data_name + '_graph.txt', dtype=np.int)<jupyter_text>Examine the Quality of the Learned GraphsCalculate the structural Hamming distance of the graphs learned by each of the algorithms<jupyter_code>shd_rai = structural_hamming_distance_cpdag(rai.graph, true_cpdag)
shd_pc = structural_hamming_distance_cpdag(pc.graph, true_cpdag)<jupyter_output><empty_output><jupyter_text>Then, calculate the Bayesian score of the true graph and the learned graph using the novel test dataset, as well as the score using the training set<jupyter_code>score_true_train = score_bdeu(true_dag, data_train, get_var_size(data_train))  # Bayesian score
score_true_test = score_bdeu(true_dag, data_test, get_var_size(data_test))

rai_dag = DAG(rai.graph.nodes_set)
rai.graph.convert_to_dag(rai_dag)
rai_score_train = score_bdeu(rai_dag, data_train, ci_test_rai.node_size)
rai_score_test = score_bdeu(rai_dag, data_test, ci_test_rai.node_size)

pc_dag = DAG(pc.graph.nodes_set)
pc.graph.convert_to_dag(pc_dag)
pc_score_train = score_bdeu(pc_dag, data_train, ci_test_pc.node_size)
pc_score_test = score_bdeu(pc_dag, data_test, ci_test_pc.node_size)<jupyter_output><empty_output><jupyter_text>Print Quality Measures of the Learned StructuresPrint Bayesian score of the graphs<jupyter_code>print('True graph BDeu: (train/test): \t{:.0f} / {:.0f}'.format(score_true_train, score_true_test))
print('RAI BDeu (train/test):  \t{:.0f} / {:.0f}'.format(rai_score_train, rai_score_test))
print('PC BDeu  (train/test):  \t{:.0f} / {:.0f}'.format(pc_score_train, pc_score_test))<jupyter_output>True graph BDeu: (train/test): 	-5060 / -48038
RAI BDeu (train/test):  	-5194 / -49219
PC BDeu  (train/test):  	-5681 / -53806<jupyter_text>Print structural hamming distance<jupyter_code>def print_shd(shd, alg_name):
    print(alg_name, '\tEdges: extra|missing', shd['edge']['extra'], shd['edge']['missing'],
          '\tOrientation: extra|missing|revered',
          shd['arrowhead']['extra'], shd['arrowhead']['missing'], shd['arrowhead']['reversed'],
          '\tTotal SHD:', shd['total'])

print_shd(shd_rai, 'RAI.')
print_shd(shd_pc, 'PC.')<jupyter_output>RAI. 	Edges: extra|missing 6 7 	Orientation: extra|missing|revered 0 3 2 	Total SHD: 18
PC. 	Edges: extra|missing 8 14 	Orientation: extra|missing|revered 0 5 0 	Total SHD: 27<jupyter_text>Plot the number of CI tests required by each algorithm<jupyter_code>num_ci_order_to_plot = 5
ci_orders = np.array(range(num_ci_order_to_plot), dtype=float)
plt.figure()
plt.bar(ci_orders-0.15, ci_test_rai.test_counter[0:num_ci_order_to_plot], width=0.3, align='center',
        color=[0, 0.25, 0.25])
plt.bar(ci_orders+0.15, ci_test_pc.test_counter[0:num_ci_order_to_plot], width=0.3, align='center',
        color=[1, 0.5, 0])
plt.title('number of CI test')
plt.xlabel('conditioning set size')
plt.legend(['RAI', 'PC'])

plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09015971148892324, 'fraction_numerical': 0.017413704276146316, 'mean_word_length': 6.358605003790751, 'pattern_counts': {'":': 0, '<': 66, '<?xml version=': 0, '>': 66, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '239836', 'license': 'Apache-2.0', 'n_tokens_mistral': 3314, 'n_tokens_neox': 3131, 'n_words': 935, 'path': 'notebooks/causal_discovery_under_causal_sufficiency.ipynb', 'repo_name': 'IntelLabs/causality-lab', 'size': '9719'}