Hi, sorry to bother you@kunwuz ,

I'm having some problems using the bic score. It always get a big number like -6951.8887543026485
rewards = local_score_BIC(data,0,[0])

The data is created by myself. I don't know the eaxt meaning of i and Pai.

Updates:

• Fix bugs (or typos) in causallearn/utils/KCI/KCI.py to conform to the original Matlab implementation in http://people.tuebingen.mpg.de/kzhang/KCI-test.zip:

  • KCI_CInd.kernel_matrix: KernelX and KernelY's empirical width depends on data_z's shape. The same applies to KernelY (Line 349, 373).
  • KCI_CInd.kernel_matrix: centering kernel matrix kzx (Line 404, 457, 462) before feeding it into KCI_V_statistic. After double checking, the optimization @MarkDana did previously (https://github.com/cmu-phil/causal-learn/pull/49) still applies.
  • KCI_CInd.kernel_matrix: change the data normalization way by adjusting the degrees of freedom correction in calculating the standard deviation.

• Add test cases where the datasets are generated by nonGaussian distributions, e.g., exponential distribution, uniform distribution and the mixture of these distributions.

• Remove print commands and add assertions on the exact p-value. The ground truth p-value is computed from the current version 26d8e36.

Resolving the concerns:

When the kernel width is set manually, there are two cases that may return a weird p-value even in the Matlab version. The details are listed in Weird_test_cases.pdf. After checking with the expert, these results are expected. More care should be taken as to what kernel width is set manually.

Comparison of the current Python version with the Matlab version on KCI-test results:

Since the Matlab version cannot support too flexible arguments, we only test the basic setting: applying Gaussian kernel to all variables X, Y and Z and setting the kernel width using empirical rules. The results for Python are as follows:

import numpy as np
import causallearn.utils.cit as cit

# np.random.seed(50)
# X = np.random.randn(30, 1)
# X_prime = np.random.randn(30, 1)
# Y = X + 0.5 * np.random.randn(30, 1)
# Z = Y + 0.5 * np.random.randn(30, 1)
# data = np.hstack((X, X_prime, Y, Z))
# np.savetxt("GaussianData.csv", data, delimiter=",")
data = np.loadtxt("GaussianData.csv", delimiter=",")

• If we use gamma approximation (approx=True), the resulting p-value is deterministic.

approx = True
for use_gp in [True, False]:
    # The argument 'use_gp' only relates the conditional KCI, so only 'pXZY' will be affected.
    cit_CIT = cit.CIT(data, 'kci', kernelX='Gaussian', kernelY='Gaussian', kernelZ='Gaussian',
                      est_width='empirical', approx=approx, use_gp=use_gp)
    pXX = cit_CIT(0, 1)
    pXZ = cit_CIT(0, 3)
    pXZY = cit_CIT(0, 3, {2})
    print('\napprox =', approx, ', use_gp =', use_gp, ':')
    print('pXX:', round(pXX, 4))
    print('pXZ:', format(pXZ, '.4e'))
    print('pXZY', round(pXZY, 4))

Output:

approx = True , use_gp = True :
pXX: 0.2662
pXZ: 5.0105e-06
pXZY 0.3585

approx = True , use_gp = False :
pXX: 0.2662
pXZ: 5.0105e-06
pXZY 0.3977

• Otherwise (approx=False), the resulting p-value is not deterministic due to the simulation of null distribution in the function null_sample_spectral(xxx) in the class KCI_UInd and KCI_CInd. (I've checked that other deterministic parts can have the same output with the Matlab version on the same input.) Therefore, we run the program 500 times and calculate the mean of p-values.

approx = False
for use_gp in [True, False]:
    np.random.seed(50)
    pXX = []
    pXZ = []
    pXZY = []
    for i in range(500):
        cit_CIT = cit.CIT(data, 'kci', kernelX='Gaussian', kernelY='Gaussian', kernelZ='Gaussian',
                  est_width='empirical', approx=approx, use_gp=use_gp)
        p01 = cit_CIT(0, 1)
        pXX.append(p01)
        p03 = cit_CIT(0, 3)
        pXZ.append(p03)
        p032 = cit_CIT(0, 3, {2})
        pXZY.append(p032)
    print('\napprox =', approx, ', use_gp =', use_gp, ':')
    print('pXX: mean = {}, var = {}'.format(round(np.mean(pXX), 4), format(np.var(pXX), '.4e')))
    print('pXZ: mean = {}, var = {}'.format(round(np.mean(pXZ), 4), format(np.var(pXZ), '.4e')))
    print('pXZY: mean = {}, var = {}'.format(round(np.mean(pXZY), 4), format(np.var(pXZY), '.4e')))

Output:

approx = False , use_gp = True :
pXX: mean = 0.246, var = 1.8853e-04
pXZ: mean = 0.0001, var = 6.6156e-08
pXZY: mean = 0.3371, var = 4.9257e-05

approx = False , use_gp = False :
pXX: mean = 0.2468, var = 1.8931e-04
pXZ: mean = 0.0001, var = 7.1376e-08
pXZY: mean = 0.3736, var = 4.5726e-05

Using the same dataset GaussianData.csv generated by Python, the Matlab version results are provided in Matlab_results.pdf.

We can see that they give the same p-value in deterministic cases and almost the same p-value in non-deterministic cases. This conclusion also applies to other datasets I tried on.

Test plan:

python -m unittest TestCIT_KCI    # should pass

TODO:

• Design more comprehensive test cases, including designing dataset of corner cases (e.g., data generated by more complicated distribution, discrete data type, larger sample size).

Hi,

I have been working with the causal-learn package for a couple of weeks now and I was wondering if there is a plan to add support for named variables?

As I understand it, at the moment data is generally assumed to take the form of a 2d np.array where each column contains a value for a particular variable X1, X2, ..., XN. I can quite easily take a csv file and convert it to this format but, in doing so, I lose the names of my columns. This relies on me having to keep track of which variable (X1, X2, ..., XN) corresponds to which named column of my csv data, which makes tasks such as visualising the output quite difficult.

I have been using the to_pyplot method to draw the discovered DAG and passing it the labels. However, this relies on the order being preserved (simply replacing X1 with the first label, X2 with the second, and so on). This has tripped me a couple of times when comparing two graphs produced from similar data where the variables do not align.

I am not sure what the solution is here, but it would be really useful if it were possible to attach actual variable names to the graph rather than doing this manually after computation. One possibility would be to support pandas dataframes for reading in data.

On my dataset (~3.6w, ~40 feats), I found it may use over 40GB memory, ( and I have not get it end successfully) Maybe it needs to use some db to reduce max memory usage?

Description section

1. add local_score_BIC_from_cov
2. set the default bic score function in GES.py with local_score_BIC_from_cov

Test plan section

python -m unittest TestGES.py # should pass

Updated functions:

• causallearn/utils/KCI/KCI.py:
  • KCI_UInd.get_kappa: reduce O(n^3) to O(n^2) based on the equation np.trace(K.dot(K)) == np.sum(K * K.T).
  • L429 elif self.kernelY == 'Polynomial': fixed a naming bug in original code. Should be kernelZ.
  • KCI_CInd.KCI_V_statistic: L479, reduced one repeated calculation based on the fact that Kzx and Kzy are usually same.
• causallearn/utils/KCI/Kernel.py:
  • Kernel.center_kernel_matrix: reduce O(n^3) to O(n^2) by plugging H = eye(n) - 1.0 / n into H.dot(K.dot(H)) and preventing the matrix multiplication.
• causallearn/search/FCMBased/lingam/hsic.py: similar trick as Kernel.center_kernel_matrix.
• causallearn/utils/cit.py: some trivial issue about specifying KCI kwargs.

Speedup gain:

• KCI_UInd (unconditional test): huge speedup (order of magnitude).
• KCI_CInd (conditional test): reduces around 30% of time.
• Overall, since conditional tests is always the bottleneck in an algorithm, this speedup is mild.
• And thus we'll need more efficient way to calculate inverse and eigens of big matrix.

Test plan:

• To ensure that this update is logically consistent with the original code (under any possible parameters):

np.random.seed(42)
from causallearn_fa79007.utils.cit import CIT as CIT_new
from causallearn_97e03ff.utils.cit import CIT as CIT_old

data = np.random.uniform(-1, 1, (100, 4))
for kernelname in ['Gaussian', 'Polynomial', 'Linear']:
    for est_width in ['empirical', 'median', 'manual']:
        for kwidth in [0.05, 0.1, 0.2]:
            for use_gp in [True, False]:
                for approx in [True, False]:
                    for polyd in [1, 2]:
                        kci_new = CIT_new(data, 'kci', kernelX=kernelname, kernelY=kernelname, kernelZ=kernelname,
                                              est_width=est_width, use_gp=use_gp, approx=approx, polyd=polyd,
                                              kwidthx=kwidth, kwidthy=kwidth, kwidthz=kwidth)
                        kci_old = CIT_old(data, 'kci', kernelX=kernelname, kernelY=kernelname, kernelZ=kernelname,
                                              est_width=est_width, use_gp=use_gp, approx=approx, polyd=polyd,
                                              kwidthx=kwidth, kwidthy=kwidth, kwidthz=kwidth)

                        # since there is randomness in null_sample_spectral, we need to fix the same seed for old and new run
                        np.random.seed(42); new_pval = kci_new(0, 1)
                        np.random.seed(42); old_pval = kci_old(0, 1)
                        assert np.isclose(old_pval, new_pval), "KCI_UIND is inconsistent after update."

                        np.random.seed(42); new_pval = kci_new(0, 1, (2, 3))
                        np.random.seed(42); old_pval = kci_old(0, 1, (2, 3))
                        assert np.isclose(old_pval, new_pval), "KCI_CIND is inconsistent after update."

• Test speedup for KCI_UInd:

for samplesize in range(1000, 20001, 1000):
    data = np.random.uniform(-1, 1, (samplesize, 2))
    kci_new = CIT_new(data, 'kci')
    kci_old = CIT_old(data, 'kci')
    tic = time.time(); kci_new(0, 1); tac = time.time(); time_new = tac - tic
    tic = time.time(); kci_old(0, 1); tac = time.time(); time_old = tac - tic
    print(f'{samplesize}:   {time_new} {time_old}')

• Test speedup for KCI_CInd:

for samplesize in range(500, 10000, 250):
    data = np.random.uniform(-1, 1, (samplesize, 4))
    kci_new = CIT_new(data, 'kci')
    kci_old = CIT_old(data, 'kci')
    tic = time.time(); kci_new(0, 1, (2, 3)); tac = time.time(); time_new = tac - tic
    tic = time.time(); kci_old(0, 1, (2, 3)); tac = time.time(); time_old = tac - tic
    print(f'{samplesize}:   {time_new} {time_old}')

After some sleuthing in the Python code and some performance testing with the help of Pablo Puig and Bryan Andrews, I've come to the conclusion that the default setting for PC in CL is maxP. That is an algorithm of mine called PCMAX. The only reference for for this algorithm currently in the literature that I know of is an arXiv tech report I put out in 2016, this one:

Ramsey, J. (2016). Improving accuracy and scalability of the pc algorithm by maximizing p-value. arXiv preprint arXiv:1610.00378.

Anyway, the performance of PC in CL is the same as the performance of PCMAX in Tetrad (and different from the performance of PC in Tetrad.)

If this is going to be made the default setting for PC for CL, perhaps this tech report should be referenced for the PC algorithm in the documentation? Otherwise no one will know what the algroithm is. I'll see if I can't sent that to a conference somewhere to get it published. I didn't realize that was being made the default here.

Update

• Add auto test to Granger Causality. The ground truth p_value_matrix_truth and adj_matrix_truth are computed from the current version b49980d

Test plan

python -m unittest tests/TestGranger.py # should pass

TODO

• Add tests for more complicated data, such as non-Gaussian data, mix of continuous and discrete data.

Description

The GIN algorithm returns the graph and the causal order. In this algorithm, the causal order and the causal graph correspond to each other, so only the causal order in the result is asserted.

Updates

• Fixed the independence test problem in GIN
• Added assertion to the previous test cases
• Removed the plot function in previous test cases for GIN
• Added tests to test GIN algorithm using the hsic independence test.
• Updated docstring
• Removed default parameter of indep_test
• Refactored the test code

Test Plan

python -m unittest tests.TestGIN # should pass

Updated files:

1. causallearn/utils/cit.py: Added exception in fisherz test
2. tests/TestCIT.py: Added a unit tests for fisherz test

Test plan section

python -m unittest tests.TestCIT # should pass

Updates

• Remove print commands and add assertions.
• Add more simulated functions and save the data.
• Fix some typos.

Test plan python -m unittest TestANM.py # should pass

Hello,

Are the edge orientation rules for the FCI implementation taken from Zhang et al (2008)? This paper introduces additional orientation rules and proves their completeness. I am wondering if the entire set of orientation rules required for completeness are implemented in this package.

1. Zhang, J. On the completeness of orientation rules for causal discovery in the presence of latent confounders and selection bias. Artificial Intelligence 172, 1873–1896 (2008).

If not, from which citation(s) are the orientation rules taken?

Thanks!

Does this repo have plans to implement the algorithm fGES[1]? fGES seems to work well for large scale problems. I wanna do some work on a large scale problem. If there is a related plan, it will help to use fGES more conveniently on the python platform, instead of calling Tetrad implemented in Java.

[1] Ramsey J, Glymour M, Sanchez-Romero R, et al. A million variables and more: the fast greedy equivalence search algorithm for learning high-dimensional graphical causal models, with an application to functional magnetic resonance images[J]. International journal of data science and analytics, 2017, 3(2): 121-129.

When will the implementation of the PNL method be completed? I saw that development has stopped since Feb 18, 2022. What are the remaining parts of the algorithm that needs to be implemented precisely? Maybe I can help

Hi,

Thank you for your great work on this package! I am testing the behaviour of the PC algorithm on simple simulated data. I found that the number of directed edges detected differs based on the ordering of the variables given to the algorithm.

I am running the following code:

import numpy as np
import matplotlib.pyplot as plt
from causallearn.search.ConstraintBased.PC import pc

def simulate_data(n_obs):
    '''
    Simulate data from the following graph

        A       B
         \     /
          v   v
            C
          /   \
         v     v
        D       E
    '''
    A = np.random.normal(size = n_obs)
    B = np.random.normal(size = n_obs)
    C = A + B + np.random.normal(size = n_obs)*0.25
    D = C + np.random.normal(size = n_obs)*0.25
    E = C + np.random.normal(size = n_obs)*0.25
    return np.stack((A,B,C,D,E), axis =1)

# generate data
n = 10000
data = simulate_data(n)

# test different permutations
permutations = [[0,1,2,3,4], [0,1,3,2,4]]

for permutation in permutations:
    graph = pc(data[:,permutation], 0.05, 'fisherz', node_names = np.array(["A","B", "C", "D", "E"])[permutation],  verbose = False)
    graph.draw_pydot_graph(labels=np.array(["A","B", "C", "D", "E"])[permutation])
    plt.show()

When the ordering of variables C and D is permuted, the PC algorithm returns the Graph with an undirected edge from C to D. However, when the ordering is unpermuted, the PC algorithm correctly directs the edge from C to D. This happens, even though the p-values in the CI tests are unchanged for the two permutations. Is this expected behaviour or can you help me fix this?

Hi, I hope to use the dag2pag function in the causal-learn lab, i.e., causallearn.utils.DAG2PAG.

However, I notice in the documentation, this function only support add latent variables.

Is there a way to add selection variables (maybe a custom change to the function) ?

Many thanks!