Causal Dynamics of Social Gaze in Primate Prefrontal-Amygdala Networks Revealed by Dynamic Bayesian Modeling.

1. from collections import deque
2. from collections.abc import Hashable
4. import networkx as nx
5. import pandas as pd
6. from sklearn.base import clone
7. from tqdm.auto import trange
9. from pgmpy import config
10. from pgmpy.base import DAG
11. from pgmpy.causal_discovery import ExpertKnowledge
12. from pgmpy.causal_discovery._base import BaseCausalDiscovery, _ScoreMixin
13. from pgmpy.structure_score import BaseStructureScore, get_scoring_method
16. class HillClimbSearch(_ScoreMixin, BaseCausalDiscovery):
17. """
18. Score-based causal discovery using hill climbing optimization.
20. This class implements the HillClimbSearch algorithm [1] for causal discovery.
21. Given a tabular dataset, the algorithm estimates the causal structure among
22. the variables in the data as a Directed Acyclic Graph (DAG). The algorithm
23. works by iteratively making local modifications to the graph structure
24. (adding, removing, or reversing edges) and keeping changes that improve
25. the score until a local maximum is reached.
27. The algorithm is a greedy local search method that:
29. 1. Starts from an initial graph (empty by default or based on provided expert knowledge).
30. 2. Evaluates all possible single-edge modifications (add, delete, reverse).
31. 3. Applies the modification with the highest score improvement.
32. 4. Repeats until no improvement can be made.
34. A tabu list is used to prevent the algorithm from immediately undoing recent
35. changes, which helps avoid getting stuck in local optima.
37. Parameters
38. ----------
39. scoring_method : str or BaseStructureScore instance, default=None
40. The score to be optimized during structure estimation. Please refer :doc:`/api/structure_score` for a list of
41. available scoring methods.
43. If ``None``, the appropriate scoring method is automatically selected based on the data type. If a string is
44. provided, the corresponding scoring method is instantiated with default parameters. To customize score-specific
45. parameters, please pass an instance of the scoring class.
47. start_dag : DAG instance, default=None
48. The starting point for the local search. By default, a completely
49. disconnected network (no edges) is used. If provided, the DAG must
50. contain exactly the same variables as in the data.
52. tabu_length : int, default=100
53. The number of recent graph modifications to store in the tabu list.
54. These modifications cannot be reversed during the search procedure.
55. This serves to enforce a wider exploration of the search space.
57. max_indegree : int or None, default=None
58. If provided, the procedure only searches among models where all nodes
59. have at most ``max_indegree`` parents. This can significantly reduce
60. the search space and computation time for large graphs.
62. expert_knowledge : ExpertKnowledge instance, default=None
63. Expert knowledge to be used with the algorithm. Expert knowledge
64. allows specification of:
66. - Required edges that must be present in the final graph
67. - Forbidden edges that cannot be present in the final graph
68. - Temporal ordering of nodes
70. return_type : str, default='pdag'
71. The type of graph to return. Options are:
73. - 'dag': Returns a directed acyclic graph (DAG).
74. - 'pdag': Returns a partially directed acyclic graph (PDAG) where edges that
75. could not be oriented are left undirected.
77. epsilon : float, default=1e-4
78. Defines the exit condition. If the improvement in score is less
79. than ``epsilon``, the algorithm terminates and returns the learned model.
81. max_iter : int, default=1e6
82. The maximum number of iterations allowed. The algorithm terminates
83. and returns the learned model when the number of iterations exceeds
84. ``max_iter``.
86. show_progress : bool, default=True
87. If True, shows a progress bar while learning the causal structure.
89. Attributes
90. ----------
91. causal_graph_ : DAG
92. The learned causal graph as a DAG at a (local) score maximum.
94. adjacency_matrix_ : pd.DataFrame
95. Adjacency matrix representation of the learned causal graph.
97. n_features_in_ : int
98. The number of features in the data used to learn the causal graph.
100. feature_names_in_ : np.ndarray
101. The feature names in the data used to learn the causal graph.
103. Examples
104. --------
105. Simulate some data to use for causal discovery:
107. >>> from pgmpy.example_models import load_model
108. >>> model = load_model("bnlearn/alarm")
109. >>> df = model.simulate(n_samples=1000, seed=42)
111. Use the HillClimbSearch algorithm to learn the causal structure from data:
113. >>> from pgmpy.causal_discovery import HillClimbSearch
114. >>> hc = HillClimbSearch(scoring_method="bic-d")
115. >>> hc.fit(df)
116. HillClimbSearch(scoring_method='bic-d')
117. >>> _ = hc.causal_graph_.edges()
119. Use expert knowledge to constrain the search:
121. >>> from pgmpy.causal_discovery import ExpertKnowledge
122. >>> expert = ExpertKnowledge(forbidden_edges=[("HISTORY", "CVP")])
123. >>> hc = HillClimbSearch(scoring_method="bic-d", expert_knowledge=expert)
124. >>> hc.fit(df) # doctest: +ELLIPSIS
125. HillClimbSearch(expert_knowledge=ExpertKnowledge(...),
126. scoring_method='bic-d')
128. References
129. ----------
130. - :footcite:t:`koller_friedman_2009`
131. """
133. def __init__(
134. self,
135. scoring_method: str | BaseStructureScore | None = None,
136. start_dag: DAG | None = None,
137. tabu_length: int = 100,
138. max_indegree: int | None = None,
139. expert_knowledge: ExpertKnowledge | None = None,
140. return_type: str = "pdag",
141. epsilon: float = 1e-4,
142. max_iter: int = int(1e6),
143. show_progress: bool = True,
144. ):
145. self.scoring_method = scoring_method
146. self.start_dag = start_dag
147. self.tabu_length = tabu_length
148. self.max_indegree = max_indegree
149. self.expert_knowledge = expert_knowledge
150. self.return_type = return_type
151. self.epsilon = epsilon
152. self.max_iter = max_iter
153. self.show_progress = show_progress
155. def _fit(self, X: pd.DataFrame):
156. """
157. The fitting procedure for the HillClimbSearch algorithm.
159. Parameters
160. ----------
161. X : pd.DataFrame or np.ndarray
162. The data to learn the causal structure from. If a numpy array is
163. passed, then the column names would be integers from 0 to n_features-1.
165. Returns
166. -------
167. self : pgmpy.causal_discovery.HillClimbSearch
168. Returns the instance with the fitted attributes.
169. """
170. self.variables_ = list(X.columns)
172. # Step 1: Initial checks and setup for arguments
173. # Step 1.1: Check score
174. score = get_scoring_method(self.scoring_method, X)
176. # Step 1.2: Check the start_dag
177. if self.start_dag is None:
178. start_dag = DAG()
179. start_dag.add_nodes_from(self.variables_)
180. elif not isinstance(self.start_dag, DAG) or not set(self.start_dag.nodes()) == set(self.variables_):
181. raise ValueError("'start_dag' should be a DAG with the same variables as the data set, or 'None'.")
182. else:
183. start_dag = self.start_dag.copy()
185. # Step 1.3: Check if expert knowledge was specified
186. if self.expert_knowledge is None:
187. expert_knowledge = ExpertKnowledge()
188. else:
189. # Clone so the fitted (`*_`) attributes land on a fresh copy, not the user's object.
190. expert_knowledge = clone(self.expert_knowledge)
192. # Step 1.3.1: Resolve the expert knowledge into its fitted (`*_`) attributes.
193. expert_knowledge.fit(X)
195. # Step 1.4: Check if required edges cause a cycle
196. start_dag.add_edges_from(expert_knowledge.required_edges_)
197. if not nx.is_directed_acyclic_graph(start_dag):
198. raise ValueError(
199. "required_edges create a cycle in start_dag. Please modify either required_edges or start_dag."
200. )
201. start_dag.remove_edges_from(expert_knowledge.forbidden_edges_)
203. # Step 1.5: Initialize max_indegree, tabu_list, and progress bar
204. max_indegree = self.max_indegree
205. if max_indegree is None:
206. max_indegree = float("inf")
208. tabu_list: deque[tuple[str, tuple[Hashable, Hashable]]] = deque(maxlen=self.tabu_length)
209. current_model = start_dag
211. if self.show_progress and config.SHOW_PROGRESS:
212. iteration = trange(int(self.max_iter))
213. else:
214. iteration = range(int(self.max_iter))
216. # Step 2: For each iteration, find the best scoring operation and
217. # do that to the current model. If no legal operation is
218. # possible, sets best_operation=None.
219. for _ in iteration:
220. best_operation, best_score_delta = max(
221. self._legal_operations_dag(
222. model=current_model,
223. scoring_method=score,
224. tabu_list=tabu_list,
225. max_indegree=max_indegree,
226. forbidden_edges=expert_knowledge.forbidden_edges_,
227. required_edges=expert_knowledge.required_edges_,
228. ),
229. key=lambda t: t[1],
230. default=(None, None),
231. )
233. if best_operation is None or best_score_delta < self.epsilon:
234. break
235. elif best_operation[0] == "+":
236. current_model.add_edge(*best_operation[1])
237. tabu_list.append(("-", best_operation[1]))
238. elif best_operation[0] == "-":
239. current_model.remove_edge(*best_operation[1])
240. tabu_list.append(("+", best_operation[1]))
241. elif best_operation[0] == "flip":
242. X_node, Y_node = best_operation[1]
243. current_model.remove_edge(X_node, Y_node)
244. current_model.add_edge(Y_node, X_node)
245. tabu_list.append(best_operation)
247. # Step 3: Store results
248. if self.return_type.lower() == "dag":
249. self.causal_graph_ = current_model
250. elif self.return_type.lower() == "pdag":
251. self.causal_graph_ = current_model.to_pdag()
252. else:
253. raise ValueError(f"return_type must be one of: dag, pdag, or cpdag. Got: {self.return_type}")
255. self.adjacency_matrix_ = self.causal_graph_.to_adjacency(
256. encoding="binary", nodelist=list(self.causal_graph_.nodes())
257. )
259. return self

HillClimbSearch.py at commit fa42ee2, under MIT · at the source