Merge branch 'py-why:main' into main

Author: kenneth-lee-ch

CODE_OF_CONDUCT.md

@@ -5,7 +5,7 @@
 We as members, contributors, and leaders pledge to make participation in our
 community a harassment-free experience for everyone, regardless of age, body
 size, visible or invisible disability, ethnicity, sex characteristics, gender
-identity and expression, level of experience, education, socio-economic status,
+identity and expression, level of experience, education, socioeconomic status,
 nationality, personal appearance, race, religion, or sexual identity
 and orientation.
 

README.md

@@ -1,6 +1,6 @@
 # causal-learn: Causal Discovery in Python
 
-Causal-learn is a python package for causal discovery that implements both classical and state-of-the-art causal discovery algorithms, which is a Python translation and extension of [Tetrad](https://github.com/cmu-phil/tetrad).
+Causal-learn ([documentation](https://causal-learn.readthedocs.io/en/latest/), [paper](https://jmlr.org/papers/volume25/23-0970/23-0970.pdf)) is a python package for causal discovery that implements both classical and state-of-the-art causal discovery algorithms, which is a Python translation and extension of [Tetrad](https://github.com/cmu-phil/tetrad).
 
 The package is actively being developed. Feedbacks (issues, suggestions, etc.) are highly encouraged.
 
@@ -76,10 +76,13 @@ Although causal-learn provides python implementations for some causal discovery
 Please cite as:
 
 ```
-@article{causallearn,
-  title={Causal-learn: Causal Discovery in Python},
-  author={Yujia Zheng and Biwei Huang and Wei Chen and Joseph Ramsey and Mingming Gong and Ruichu Cai and Shohei Shimizu and Peter Spirtes and Kun Zhang},
-  journal={arXiv preprint arXiv:2307.16405},
-  year={2023}
+@article{zheng2024causal,
+  title={Causal-learn: Causal discovery in python},
+  author={Zheng, Yujia and Huang, Biwei and Chen, Wei and Ramsey, Joseph and Gong, Mingming and Cai, Ruichu and Shimizu, Shohei and Spirtes, Peter and Zhang, Kun},
+  journal={Journal of Machine Learning Research},
+  volume={25},
+  number={60},
+  pages={1--8},
+  year={2024}
 }
 ```

causallearn/score/LocalScoreFunction.py

@@ -34,9 +34,9 @@ def local_score_BIC(Data: ndarray, i: int, PAi: List[int], parameters=None) -> f
     if len(PAi) == 0:
         return n * np.log(cov[i, i])
 
-    yX = np.mat(cov[np.ix_([i], PAi)])
-    XX = np.mat(cov[np.ix_(PAi, PAi)])
-    H = np.log(cov[i, i] - yX * np.linalg.inv(XX) * yX.T)
+    yX = cov[np.ix_([i], PAi)]
+    XX = cov[np.ix_(PAi, PAi)]
+    H = np.log(cov[i, i] - yX @ np.linalg.inv(XX) @ yX.T)
 
     return n * H + np.log(n) * len(PAi) * lambda_value
 
@@ -68,9 +68,9 @@ def local_score_BIC_from_cov(
     if len(PAi) == 0:
         return n * np.log(cov[i, i])
 
-    yX = np.mat(cov[np.ix_([i], PAi)])
-    XX = np.mat(cov[np.ix_(PAi, PAi)])
-    H = np.log(cov[i, i] - yX * np.linalg.inv(XX) * yX.T)
+    yX = cov[np.ix_([i], PAi)]
+    XX = cov[np.ix_(PAi, PAi)]
+    H = np.log(cov[i, i] - yX @ np.linalg.inv(XX) @ yX.T)
 
     return n * H + np.log(n) * len(PAi) * lambda_value
 

causallearn/score/LocalScoreFunctionClass.py

@@ -46,3 +46,9 @@ def score(self, i: int, PAi: List[int]) -> float:
                 self.score_cache[i][hash_key] = self.local_score_fun(self.data, i, PAi, self.parameters)
 
         return self.score_cache[i][hash_key]
+
+    def score_nocache(self, i: int, PAi: List[int]) -> float:
+        if self.local_score_fun == local_score_BIC_from_cov:
+            return self.local_score_fun((self.cov, self.n), i, PAi, self.parameters)
+        else:
+            return self.local_score_fun(self.data, i, PAi, self.parameters)

causallearn/search/ConstraintBased/FCI.py

@@ -320,8 +320,9 @@ def rulesR1R2cycle(graph: Graph, bk: BackgroundKnowledge | None, changeFlag: boo
 
 def isNoncollider(graph: Graph, sep_sets: Dict[Tuple[int, int], Set[int]], node_i: Node, node_j: Node,
                   node_k: Node) -> bool:
-    sep_set = sep_sets[(graph.get_node_map()[node_i], graph.get_node_map()[node_k])]
-    return sep_set is not None and sep_set.__contains__(graph.get_node_map()[node_j])
+    node_map = graph.get_node_map()
+    sep_set = sep_sets.get((node_map[node_i], node_map[node_k]))
+    return sep_set is not None and sep_set.__contains__(node_map[node_j])
 
 
 def ruleR3(graph: Graph, sep_sets: Dict[Tuple[int, int], Set[int]], bk: BackgroundKnowledge | None, changeFlag: bool,
@@ -827,7 +828,8 @@ def _contains_all(set_a: Set[Node], set_b: Set[Node]):
 
 
 def fci(dataset: ndarray, independence_test_method: str=fisherz, alpha: float = 0.05, depth: int = -1,
-        max_path_length: int = -1, verbose: bool = False, background_knowledge: BackgroundKnowledge | None = None, show_progress: bool = True,
+        max_path_length: int = -1, verbose: bool = False, background_knowledge: BackgroundKnowledge | None = None, 
+        show_progress: bool = True, node_names = None,
         **kwargs) -> Tuple[Graph, List[Edge]]:
     """
     Perform Fast Causal Inference (FCI) algorithm for causal discovery
@@ -882,8 +884,10 @@ def fci(dataset: ndarray, independence_test_method: str=fisherz, alpha: float =
 
 
     nodes = []
+    if node_names is None:
+        node_names = [f"X{i + 1}" for i in range(dataset.shape[1])]
     for i in range(dataset.shape[1]):
-        node = GraphNode(f"X{i + 1}")
+        node = GraphNode(node_names[i])
         node.add_attribute("id", i)
         nodes.append(node)
 

causallearn/search/FCMBased/lingam/var_lingam.py

@@ -319,10 +319,14 @@ def _pruning(self, X, B_taus, causal_order):
             obj = np.zeros((len(blocks)))
             exp = np.zeros(
                 (len(blocks), causal_order_no + n_features * self._lags))
-            for j, block in enumerate(blocks):
-                obj[j] = block[0][i]
-                exp[j:] = np.concatenate(
-                    [block[0][ancestor_indexes].flatten(), block[1:][:].flatten()], axis=0)
+            
+            # Create 3D array to hold flattened ancestor indices for each block
+            ancestor_indexes_flat = blocks[:, 0, ancestor_indexes].reshape(len(blocks), -1)
+            # Fill obj using advanced indexing
+            obj[:] = blocks[:, 0, i]
+            # Fill exp using advanced indexing
+            exp[:, :causal_order_no] = ancestor_indexes_flat
+            exp[:, causal_order_no:] = blocks[:, 1:].reshape(len(blocks), -1)
 
             # adaptive lasso
             gamma = 1.0

causallearn/search/PermutationBased/BOSS.py

@@ -0,0 +1,207 @@
+import random
+import sys
+import time
+import warnings
+from typing import Any, Dict, List, Optional
+
+import numpy as np
+from causallearn.graph.GeneralGraph import GeneralGraph
+from causallearn.graph.GraphNode import GraphNode
+from causallearn.score.LocalScoreFunction import (
+    local_score_BDeu,
+    local_score_BIC,
+    local_score_BIC_from_cov,
+    local_score_cv_general,
+    local_score_cv_multi,
+    local_score_marginal_general,
+    local_score_marginal_multi,
+)
+from causallearn.search.PermutationBased.gst import GST;
+from causallearn.score.LocalScoreFunctionClass import LocalScoreClass
+from causallearn.utils.DAG2CPDAG import dag2cpdag
+
+
+def boss(
+    X: np.ndarray,
+    score_func: str = "local_score_BIC_from_cov",
+    parameters: Optional[Dict[str, Any]] = None,
+    verbose: Optional[bool] = True,
+    node_names: Optional[List[str]] = None,
+) -> GeneralGraph:
+    """
+    Perform a best order score search (BOSS) algorithm
+
+    Parameters
+    ----------
+    X : data set (numpy ndarray), shape (n_samples, n_features). The input data, where n_samples is the number of samples and n_features is the number of features.
+    score_func : the string name of score function. (str(one of 'local_score_CV_general', 'local_score_marginal_general',
+                    'local_score_CV_multi', 'local_score_marginal_multi', 'local_score_BIC', 'local_score_BIC_from_cov', 'local_score_BDeu')).
+    parameters : when using CV likelihood,
+                  parameters['kfold']: k-fold cross validation
+                  parameters['lambda']: regularization parameter
+                  parameters['dlabel']: for variables with multi-dimensions,
+                               indicate which dimensions belong to the i-th variable.
+    verbose : whether to print the time cost and verbose output of the algorithm.
+
+    Returns
+    -------
+    G : learned causal graph, where G.graph[j,i] = 1 and G.graph[i,j] = -1 indicates i --> j, G.graph[i,j] = G.graph[j,i] = -1 indicates i --- j.
+    """
+
+    X = X.copy()
+    n, p = X.shape
+    if n < p:
+        warnings.warn("The number of features is much larger than the sample size!")
+
+    if score_func == "local_score_CV_general":  
+        # % k-fold negative cross validated likelihood based on regression in RKHS
+        if parameters is None:
+            parameters = {
+                "kfold": 10,  # 10 fold cross validation
+                "lambda": 0.01,
+            }  # regularization parameter
+        localScoreClass = LocalScoreClass(
+            data=X, local_score_fun=local_score_cv_general, parameters=parameters
+        )
+    elif score_func == "local_score_marginal_general":  
+        # negative marginal likelihood based on regression in RKHS
+        parameters = {}
+        localScoreClass = LocalScoreClass(
+            data=X, local_score_fun=local_score_marginal_general, parameters=parameters
+        )
+    elif score_func == "local_score_CV_multi":  
+        # k-fold negative cross validated likelihood based on regression in RKHS
+        # for data with multi-variate dimensions
+        if parameters is None:
+            parameters = {
+                "kfold": 10,
+                "lambda": 0.01,
+                "dlabel": {},
+            }  # regularization parameter
+            for i in range(X.shape[1]):
+                parameters["dlabel"]["{}".format(i)] = i
+        localScoreClass = LocalScoreClass(
+            data=X, local_score_fun=local_score_cv_multi, parameters=parameters
+        )
+    elif score_func == "local_score_marginal_multi":  
+        # negative marginal likelihood based on regression in RKHS
+        # for data with multi-variate dimensions
+        if parameters is None:
+            parameters = {"dlabel": {}}
+            for i in range(X.shape[1]):
+                parameters["dlabel"]["{}".format(i)] = i
+        localScoreClass = LocalScoreClass(
+            data=X, local_score_fun=local_score_marginal_multi, parameters=parameters
+        )
+    elif score_func == "local_score_BIC":  
+        # SEM BIC score
+        warnings.warn("Please use 'local_score_BIC_from_cov' instead")
+        if parameters is None:
+            parameters = {"lambda_value": 2}
+        localScoreClass = LocalScoreClass(
+            data=X, local_score_fun=local_score_BIC, parameters=parameters
+        )
+    elif score_func == "local_score_BIC_from_cov":  
+        # SEM BIC score
+        if parameters is None:
+            parameters = {"lambda_value": 2}
+        localScoreClass = LocalScoreClass(
+            data=X, local_score_fun=local_score_BIC_from_cov, parameters=parameters
+        )
+    elif score_func == "local_score_BDeu":  
+        # BDeu score
+        localScoreClass = LocalScoreClass(
+            data=X, local_score_fun=local_score_BDeu, parameters=None
+        )
+    else:
+        raise Exception("Unknown function!")
+    
+    score = localScoreClass
+    gsts = [GST(i, score) for i in range(p)]
+
+    node_names = [("X%d" % (i + 1)) for i in range(p)] if node_names is None else node_names
+    nodes = []
+
+    for name in node_names:
+        node = GraphNode(name)
+        nodes.append(node)
+
+    G = GeneralGraph(nodes)
+
+    runtime = time.perf_counter()
+        
+    order = [v for v in range(p)]
+
+    gsts = [GST(v, score) for v in order]
+    parents = {v: [] for v in order}
+    
+    variables = [v for v in order]
+    while True:
+        improved = False
+        random.shuffle(variables)
+        if verbose:
+            for i, v in enumerate(order):
+                parents[v].clear()
+                gsts[v].trace(order[:i], parents[v])
+            sys.stdout.write("\rBOSS edge count: %i    " % np.sum([len(parents[v]) for v in range(p)]))
+            sys.stdout.flush()
+
+        for v in variables:
+            improved |= better_mutation(v, order, gsts)
+        if not improved: break
+
+    for i, v in enumerate(order):
+        parents[v].clear()
+        gsts[v].trace(order[:i], parents[v])
+
+    runtime = time.perf_counter() - runtime
+    
+    if verbose:
+        sys.stdout.write("\nBOSS completed in: %.2fs \n" % runtime)
+        sys.stdout.flush()
+
+    for y in range(p):
+        for x in parents[y]:
+            G.add_directed_edge(nodes[x], nodes[y])
+
+    G = dag2cpdag(G)
+
+    return G
+
+
+def reversed_enumerate(iter, j):
+    for w in reversed(iter):
+        yield j, w
+        j -= 1
+
+
+def better_mutation(v, order, gsts):
+    i = order.index(v)
+    p = len(order)
+    scores = np.zeros(p + 1)
+
+    prefix = []
+    score = 0
+    for j, w in enumerate(order):
+        scores[j] = gsts[v].trace(prefix) + score
+        if v != w:
+            score += gsts[w].trace(prefix)
+            prefix.append(w)
+
+    scores[p] = gsts[v].trace(prefix) + score
+    best = p
+
+    prefix.append(v)
+    score = 0
+    for j, w in reversed_enumerate(order, p - 1):
+        if v != w:
+            prefix.remove(w)
+            score += gsts[w].trace(prefix)
+        scores[j] += score
+        if scores[j] > scores[best]: best = j
+        
+    if scores[i] + 1e-6 > scores[best]: return False
+    order.remove(v)
+    order.insert(best - int(best > i), v)
+
+    return True