因果推理算法及工具应用

from notears import linear, utils
import numpy as np
from tqdm import tqdm
from collections import defaultdict
import os


def main():
    utils.set_random_seed(123)
    # 生成图的数量
    num_graph = 2
    # 每个图对应生成数据集的份数
    num_data_per_graph = 2
    # 生成数据的保存地址
    save_data_path = r"data"
    # n-生成的数据量，多少条数据；d-生成图/生成数据的维度，几个变量；s0-生成图的边数；
    # graph_type-生成图的类型（ER, SF, BP）；sem_type-生成数据的分布类型（线性：gauss, exp, gumbel, uniform, logistic, poisson/非线性：mlp, mim, gp, gp-add）
    n, d, s0, graph_type, sem_type = 200, 7, 14, 'ER', 'gauss'
    n1, d1, s1, graph_type1, sem_type1 = 30, 3, 2, 'SF', 'mlp'

    # 线性数据例子
    # w_ranges-参数的数据范围
    w_ranges = ((-2.0, -0.5), (0.5, 2.0))
    # noise_scale-加性噪声的范围
    noise_scale = np.ones(d)
    # 生成数据存放文件夹的名称
    expt_name = 'equal_var'
    run_expt(save_data_path, num_graph, num_data_per_graph, n, d, s0, graph_type, sem_type, w_ranges, noise_scale, expt_name)

    # 非线性数据例子
    w_ranges = ((-2.0, -1.1), (1.1, 2.0))
    noise_scale1 = [1., 0.15, 0.2]
    expt_name = 'large_a'
    #run_expt_(save_data_path, num_graph, num_data_per_graph, n1, d1, s1, graph_type1, sem_type1, w_ranges, noise_scale1, expt_name)

# 线性
def run_expt(save_data_path, num_graph, num_data_per_graph, n, d, s0, graph_type, sem_type, w_ranges, noise_scale, expt_name):
    # os.mkdir(expt_name)
    # os.chmod(expt_name, 0o777)
    # perf = defaultdict(list)
    noise_scale  =np.array(noise_scale)
    for ii in tqdm(range(num_graph)):
        # B_true-生成因果图的0-1邻接矩阵
        B_true = utils.simulate_dag(d, s0, graph_type)
        W_true = utils.simulate_parameter(B_true, w_ranges=w_ranges)
        B_true_fn = os.path.join(save_data_path, expt_name, f'graph{ii:05}_W_true.csv')
        np.savetxt(B_true_fn, B_true, delimiter=',')
        for jj in range(num_data_per_graph):
            # X = utils.simulate_linear_sem(W_true, n, sem_type, noise_scale=noise_scale)
            X = utils.simulate_linear_sem(W_true, n, sem_type, noise_scale=noise_scale)
            X_fn = os.path.join(save_data_path, expt_name, f'graph{ii:05}_data{jj:05}_X.csv')
            np.savetxt(X_fn, X, delimiter=',')

# 非线性
def run_expt_(save_data_path, num_graph, num_data_per_graph, n, d, s0, graph_type, sem_type, noise_scale, expt_name):
    noise_scale = np.array(noise_scale)
    for ii in tqdm(range(num_graph)):
        B_true = utils.simulate_dag(d, s0, graph_type)
        B_true_fn = os.path.join(save_data_path, expt_name, f'graph{ii:05}_W_true.csv')
        np.savetxt(B_true_fn, B_true, delimiter=',')
        for jj in range(num_data_per_graph):
            X = utils.simulate_nonlinear_sem(B_true, n, sem_type, noise_scale=noise_scale)
            X_fn = os.path.join(save_data_path, expt_name, f'graph{ii:05}_data{jj:05}_X.csv')
            np.savetxt(X_fn, X, delimiter=',')


if __name__ == '__main__':
    main()

from notears import linear, nonlinear, utils
import numpy as np

with open(r'D:\Desktop\TBFILE\学习\群体智能\群体智能第二次实验\notears-master\data\equal_var\graph00001_W_true.csv') as f:
    line = f.readline()
    data_array = []
    while line:
        num = list(map(float, line.split(',')))
        data_array.append(num)
        line = f.readline()
    B_true = np.array(data_array)

with open(r'D:\Desktop\TBFILE\学习\群体智能\群体智能第二次实验\notears-master\data\equal_var\graph00001_data00001_X.csv') as f:
    line = f.readline()
    data_array = []
    while line:
        num = list(map(float, line.split(',')))
        data_array.append(num)
        line = f.readline()
    X = np.array(data_array)

param_dist = {
              # "max_iter":[100, 200, 300, 400, 500],
              "w_threshold":[0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95],
              # "lambda0":[0.001,  0.003, 0.004, 0.005,  0.007, 0.009, 0.011, 0.013, 0.015, 0.017, 0.019, 0.021, 0.023, 0.025],
              "lambda1":[0.01,  0.02, 0.03, 0.04,  0.05, 0.06, 0.07, 0.08, 0.09, 0.1],
              }

# MI = param_dist.get("max_iter")
WT = param_dist.get("w_threshold")
L1 = param_dist.get("lambda1")

result = []
for i in range(len(L1)):
    for j in range(len(WT)):
        #notears
        lambda1 = L1[i]
        w_threshold = WT[j]
        W_notears = linear.notears_linear(X, lambda1=lambda1, loss_type='l2',max_iter=300,
                                          w_threshold=w_threshold)

        assert utils.is_dag(W_notears)
        #eval
        B_notears = (W_notears !=0)
        acc = utils.count_accuracy(B_true, B_notears)
        print("----------------------------------")
        print("lambda1:",L1[i],"w_threshold:",WT[j])
        ret = [L1[i],WT[j]]

        v_met = list(acc.values())
        print("fdr|tpr|shd|nnz")
        print(v_met)
        ret += v_met
        result.append(ret)

import csv

# 创建csv writer对象
with open(r'data\equal_var\res_save\r_l11.csv', 'w', newline='') as file1:
    writer1 = csv.writer(file1)

    # 写入数据
    for row in result:
        writer1.writerow(row)

# 数据
data = [["lambda1:","w_threshold:"] + list(acc.keys())]+result

# 创建csv writer对象
with open(r'data\equal_var\res_save\d_l11.csv', 'w', newline='') as file2:
    writer2 = csv.writer(file2)

    # 写入数据
    for row in data:
        writer2.writerow(row)

import torch
import torch.nn as nn
import numpy as np
import csv
from notears.locally_connected import LocallyConnected
from notears.lbfgsb_scipy import LBFGSBScipy
from notears.trace_expm import trace_expm
import notears.utils as ut
from notears import linear, nonlinear, utils


# Define NotearsMLP and NotearsSobolev classes here (as provided)

def squared_loss(output, target):
    n = target.shape[0]
    loss = 0.5 / n * torch.sum((output - target) ** 2)
    return loss


def dual_ascent_step(model, X, lambda1, lambda2, rho, alpha, h, rho_max):
    h_new = None
    optimizer = LBFGSBScipy(model.parameters())
    X_torch = torch.from_numpy(X)
    while rho < rho_max:
        def closure():
            optimizer.zero_grad()
            X_hat = model(X_torch)
            loss = squared_loss(X_hat, X_torch)
            h_val = model.h_func()
            penalty = 0.5 * rho * h_val * h_val + alpha * h_val
            l2_reg = 0.5 * lambda2 * model.l2_reg()
            l1_reg = lambda1 * model.fc1_l1_reg()
            primal_obj = loss + penalty + l2_reg + l1_reg
            primal_obj.backward()
            return primal_obj

        optimizer.step(closure)
        with torch.no_grad():
            h_new = model.h_func().item()
        if h_new > 0.25 * h:
            rho *= 10
        else:
            break
    alpha += rho * h_new
    return rho, alpha, h_new


def notears_nonlinear(model: nn.Module,
                      X: np.ndarray,
                      lambda1: float = 0.,
                      lambda2: float = 0.,
                      max_iter: int = 100,
                      h_tol: float = 1e-8,
                      rho_max: float = 1e+16,
                      w_threshold: float = 0.3):
    rho, alpha, h = 1.0, 0.0, np.inf
    for _ in range(max_iter):
        rho, alpha, h = dual_ascent_step(model, X, lambda1, lambda2,
                                         rho, alpha, h, rho_max)
        if h <= h_tol or rho >= rho_max:
            break
    W_est = model.fc1_to_adj()
    W_est[np.abs(W_est) < w_threshold] = 0
    return W_est


def main():
    torch.set_default_dtype(torch.double)
    np.set_printoptions(precision=3)
    ut.set_random_seed(123)

    # Load observational data from CSV file
    X = np.genfromtxt(
        r'D:\Desktop\TBFILE\学习\群体智能\群体智能第二次实验\exp_nonlinear_3\graph00001_data00000_X.csv',
        delimiter=',', skip_header=1)
    # Load ground truth DAG from CSV file
    B_true = np.genfromtxt(
        r'D:\Desktop\TBFILE\学习\群体智能\群体智能第二次实验\exp_nonlinear_3\graph00001_W_true.csv',
        delimiter=',')

    d = X.shape[1]

    # Hyperparameters for grid search
    param_dist = {
        "w_threshold": [0.05,0.10,0.15,0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5],
        "lambda1": [0.01, 0.02, 0.03,  0.04, 0.05],
        "lambda2": [0.01, 0.02, 0.03, 0.04, 0.05]
    }

    WT = param_dist["w_threshold"]
    L1 = param_dist["lambda1"]
    L2 = param_dist["lambda2"]

    result = []
    for lambda1 in L1:
        for lambda2 in L2:
            for w_threshold in WT:
                model = nonlinear.NotearsMLP(dims=[d, 10, 1], bias=True)
                W_est = notears_nonlinear(model, X, lambda1=lambda1, lambda2=lambda2, w_threshold=w_threshold)

                # Ensure W_est is a DAG
                if not ut.is_dag(W_est):
                    continue  # Skip this combination if it's not a DAG

                # Save results for each combination
                ret = [lambda1, lambda2, w_threshold]
                ret += list(ut.count_accuracy(B_true, W_est != 0).values())
                result.append(ret)

                # Print results for debugging
                print(f"lambda1: {lambda1}, lambda2: {lambda2}, w_threshold: {w_threshold}")
                print("fdr|tpr|shd|nnz")
                print(list(ut.count_accuracy(B_true, W_est != 0).values()))

    # Save results to CSV files
    headers = ["lambda1", "lambda2", "w_threshold"] + list(ut.count_accuracy(B_true, W_est != 0).keys())
    data = [headers] + result

    with open(
            r'D:\Desktop\TBFILE\学习\群体智能\群体智能第二次实验\exp_nonlinear_3\grid_search_results1.csv',
            'w', newline='') as file:
        writer = csv.writer(file)
        for row in data:
            writer.writerow(row)


if __name__ == '__main__':
    main()
    print('done')

from castle.common import GraphDAG
from castle.metrics import MetricsDAG
from castle.algorithms import GraNDAG
import numpy as np
import pandas as pd
import itertools

# Set random seed for reproducibility
np.random.seed(42)

# Load observational data from CSV file
X = np.genfromtxt(
    r'D:\Desktop\TBFILE\学习\群体智能\群体智能第二次实验\exp_nonlinear_1\graph00001_data00000_X.csv',
    delimiter=',', skip_header=1)

# Load ground truth DAG from CSV file
B_true = np.genfromtxt(
    r'D:\Desktop\TBFILE\学习\群体智能\群体智能第二次实验\exp_nonlinear_1\graph00001_W_true.csv',
    delimiter=',')

# Normalize data
normalize = True
if normalize:
    X_mean = np.mean(X, axis=0, keepdims=True)
    X_std = np.std(X, axis=0, keepdims=True)
    X = (X - X_mean) / X_std

# Check for NaN values in the data
if np.isnan(X).any():
    raise ValueError("Input data contains NaN values. Please check the data preprocessing steps.")

# Define hyperparameter grid
param_grid = {
    'model_name': ['NonLinGauss'],
    'nonlinear': ['leaky-relu'],
    'optimizer': ['sgd'],
    'norm_prod': ['paths'],
    'device_type': ['cpu'],
    'random_seed': [42],
    'normalize': [normalize],
    'input_dim': [X.shape[1]],
    'hidden_num': [2, 3,4],
    'hidden_dim': [10, 20, 30],
    'batch_size': [64,128],
    'lr': [0.0001],
    'iterations': [1000]
}

# Create a list to store results
results = []

# Perform grid search
keys, values = zip(*param_grid.items())
for v in itertools.product(*values):
    params = dict(zip(keys, v))

    try:
        # Create and train GraNDAG instance
        gnd = GraNDAG(**params)
        gnd.learn(data=X)

        # Evaluate the learned DAG
        mm = MetricsDAG(gnd.causal_matrix, B_true)

        # Store results
        result = params.copy()
        result.update(mm.metrics)
        results.append(result)

        # Print current result
        print(f"Evaluated params: {params}")
        print(f"Metrics: {mm.metrics}\n")

    except ValueError as e:
        print(f"Error with params: {params}")
        print(e)
        continue

# Save results to CSV
results_df = pd.DataFrame(results)
results_df.to_csv(r'D:\Desktop\TBFILE\学习\群体智能\群体智能第二次实验\exp_nonlinear_1\pc_results.csv',
                  index=False)

print("Grid search completed. Results saved to CSV.")

import pandas as pd
import numpy as np
from castle.algorithms import PC
from castle.metrics import MetricsDAG
from bayes_opt import BayesianOptimization
import logging

# Set up logging
logging.basicConfig(level=logging.INFO)
logger = logging.getLogger(__name__)

# 读取CSV数据文件
data_path = r'D:\Desktop\TBFILE\学习\群体智能\群体智能第二次实验\exp_nonlinear_3\graph00001_data00000_X.csv'
data = pd.read_csv(data_path)

# 读取真实的因果图
true_graph_path = r'D:\Desktop\TBFILE\学习\群体智能\群体智能第二次实验\exp_nonlinear_3\graph00001_W_true.csv'
true_dag = pd.read_csv(true_graph_path, header=None).values

# 将数据转换为numpy数组
X = data.values

# 定义目标函数
def objective(alpha):
    try:
        # 设置新的alpha值并初始化PC算法
        pc = PC(variant='stable', alpha=alpha)

        # 学习因果结构
        pc.learn(X)

        # 评估学习到的因果矩阵
        met = MetricsDAG(pc.causal_matrix, true_dag)

        # 记录SHD值
        shd = met.metrics['shd']
        
        logger.info(f'Alpha: {alpha}, SHD: {shd}')

        # 返回SHD的负值，因为贝叶斯优化是一个最大化器，我们希望最小化SHD
        return -shd
    except Exception as e:
        logger.error(f"Error with alpha: {alpha}")
        logger.error(e)
        return -float('inf')

# 设置贝叶斯优化的参数范围
pbounds = {'alpha': (0.01, 0.1)}

# 初始化贝叶斯优化器
optimizer = BayesianOptimization(
    f=objective,
    pbounds=pbounds,
    random_state=42
)

# 执行优化过程
optimizer.maximize(n_iter=20)

# 获取最佳参数
best_alpha = optimizer.max['params']['alpha']
best_shd = -optimizer.max['target']

print(f'\nBest alpha: {best_alpha}, with SHD: {best_shd}')