diff --git a/examples/gcil/README.md b/examples/gcil/README.md
new file mode 100644
index 00000000..ac7ca10b
--- /dev/null
+++ b/examples/gcil/README.md
@@ -0,0 +1,35 @@
+# Graph Contrastive Invariant Learning from the Causal Perspective (GCIL)
+
+- Paper link: [https://arxiv.org/pdf/2401.12564v2](https://arxiv.org/pdf/2401.12564v2)
+- Author's code repo: [https://github.com/BUPT-GAMMA/GCIL](https://github.com/BUPT-GAMMA/GCIL). Note that the original code is 
+  implemented with Tensorflow for the paper. 
+
+# Dataset Statics
+
+| Dataset | # Nodes | # Edges | # Classes |
+| ------- | ------- | ------- | --------- |
+| Cora    | 2,708   | 10,556  | 7         |
+| Pubmed  | 19,717  | 88,651  | 3         |
+
+Refer to [Planetoid](https://gammagl.readthedocs.io/en/latest/api/gammagl.datasets.html#gammagl.datasets.Planetoid).
+
+Results
+-------
+
+```bash
+# available dataset: "cora", "pubmed"
+TL_BACKEND="torch" python gcil_trainer.py cora
+TL_BACKEND="torch" python gcil_trainer.py pubmed 
+```
+
+Ma-F1:
+| Dataset | Paper    | Our(th)    |
+| ------- | -------- | ---------- |
+| cora    | 83.8±0.5 | 45.19±0.22 |
+| pubmed  | 81.5±0.5 | 46.30±0.02 |
+
+Mi-F1
+| Dataset | Paper    | Our(th)    |
+| ------- | -------- | ---------- |
+| cora    | 84.4±0.7 | 49.71±0.22 |
+| pubmed  | 81.6±0.7 | 53.77±0.01 |
\ No newline at end of file
diff --git a/examples/gcil/aug.py b/examples/gcil/aug.py
new file mode 100644
index 00000000..e52a1be5
--- /dev/null
+++ b/examples/gcil/aug.py
@@ -0,0 +1,50 @@
+import torch as th
+import numpy as np
+from gammagl.data import Graph
+
+
+def random_aug(graph, attr, diag_attr, x, feat_drop_rate, edge_mask_rate, device):
+    n_node = graph.num_nodes if hasattr(graph, 'num_nodes') else graph.x.shape[0]
+
+    edge_mask = mask_edge(graph, edge_mask_rate)
+
+    feat = drop_feature(x, feat_drop_rate)
+
+    edge_index = graph.edge_index[:, edge_mask].long()  
+
+    edge_weight = attr[edge_mask] if attr is not None else None
+
+    if isinstance(attr, np.ndarray):
+        attr_cpu = attr.cpu().numpy() if attr.is_cuda else attr.numpy()
+        diag_attr_cpu = diag_attr.cpu().numpy() if diag_attr.is_cuda else diag_attr.numpy()
+        attr = np.concatenate([attr_cpu[edge_mask], diag_attr_cpu], axis=0)
+
+    new_graph = Graph(x=feat, edge_index=edge_index)
+
+    new_graph.x = new_graph.x.to(device)
+    new_graph.edge_index = new_graph.edge_index.to(device)
+
+    if edge_weight is not None:
+        edge_weight = th.tensor(edge_weight, dtype=th.float32).to(device)  # 确保 edge_weight 是 FloatTensor
+
+    return new_graph, edge_weight, feat
+
+
+
+def drop_feature(x, drop_prob):
+    drop_mask = th.empty(
+        (x.size(1),),
+        dtype=th.float32,
+        device=x.device).uniform_(0, 1) < drop_prob
+    x = x.clone()
+    x[:, drop_mask] = 0
+
+    return x
+
+
+def mask_edge(graph, edge_mask_rate):
+    E = graph.edge_index.shape[1]
+
+    mask = np.random.rand(E) > edge_mask_rate
+
+    return mask
\ No newline at end of file
diff --git a/examples/gcil/dataset/cora/0.01_1_0.npz b/examples/gcil/dataset/cora/0.01_1_0.npz
new file mode 100644
index 00000000..812b6ca1
Binary files /dev/null and b/examples/gcil/dataset/cora/0.01_1_0.npz differ
diff --git a/examples/gcil/dataset/cora/0.01_1_1.npz b/examples/gcil/dataset/cora/0.01_1_1.npz
new file mode 100644
index 00000000..2ca20a3a
Binary files /dev/null and b/examples/gcil/dataset/cora/0.01_1_1.npz differ
diff --git a/examples/gcil/dataset/cora/0.01_1_2.npz b/examples/gcil/dataset/cora/0.01_1_2.npz
new file mode 100644
index 00000000..b748e8f3
Binary files /dev/null and b/examples/gcil/dataset/cora/0.01_1_2.npz differ
diff --git a/examples/gcil/dataset/cora/0.01_1_3.npz b/examples/gcil/dataset/cora/0.01_1_3.npz
new file mode 100644
index 00000000..8ea99b6f
Binary files /dev/null and b/examples/gcil/dataset/cora/0.01_1_3.npz differ
diff --git a/examples/gcil/dataset/cora/0.01_1_4.npz b/examples/gcil/dataset/cora/0.01_1_4.npz
new file mode 100644
index 00000000..3ae63fd3
Binary files /dev/null and b/examples/gcil/dataset/cora/0.01_1_4.npz differ
diff --git a/examples/gcil/dataset/cora/0.01_1_5.npz b/examples/gcil/dataset/cora/0.01_1_5.npz
new file mode 100644
index 00000000..0c0ebe1f
Binary files /dev/null and b/examples/gcil/dataset/cora/0.01_1_5.npz differ
diff --git a/examples/gcil/dataset/cora/0.01_1_6.npz b/examples/gcil/dataset/cora/0.01_1_6.npz
new file mode 100644
index 00000000..2b0ce480
Binary files /dev/null and b/examples/gcil/dataset/cora/0.01_1_6.npz differ
diff --git a/examples/gcil/dataset/pubmed/0.01_1_0.npz b/examples/gcil/dataset/pubmed/0.01_1_0.npz
new file mode 100644
index 00000000..c83cf7b3
Binary files /dev/null and b/examples/gcil/dataset/pubmed/0.01_1_0.npz differ
diff --git a/examples/gcil/dataset/pubmed/0.01_1_1.npz b/examples/gcil/dataset/pubmed/0.01_1_1.npz
new file mode 100644
index 00000000..40040f3b
Binary files /dev/null and b/examples/gcil/dataset/pubmed/0.01_1_1.npz differ
diff --git a/examples/gcil/dataset/pubmed/0.01_1_2.npz b/examples/gcil/dataset/pubmed/0.01_1_2.npz
new file mode 100644
index 00000000..07af3567
Binary files /dev/null and b/examples/gcil/dataset/pubmed/0.01_1_2.npz differ
diff --git a/examples/gcil/dataset/pubmed/0.01_1_3.npz b/examples/gcil/dataset/pubmed/0.01_1_3.npz
new file mode 100644
index 00000000..4b0414ca
Binary files /dev/null and b/examples/gcil/dataset/pubmed/0.01_1_3.npz differ
diff --git a/examples/gcil/dataset/pubmed/0.01_1_4.npz b/examples/gcil/dataset/pubmed/0.01_1_4.npz
new file mode 100644
index 00000000..2420d266
Binary files /dev/null and b/examples/gcil/dataset/pubmed/0.01_1_4.npz differ
diff --git a/examples/gcil/dataset/pubmed/0.01_1_5.npz b/examples/gcil/dataset/pubmed/0.01_1_5.npz
new file mode 100644
index 00000000..04d95c45
Binary files /dev/null and b/examples/gcil/dataset/pubmed/0.01_1_5.npz differ
diff --git a/examples/gcil/dataset/pubmed/0.01_1_6.npz b/examples/gcil/dataset/pubmed/0.01_1_6.npz
new file mode 100644
index 00000000..5fa79be8
Binary files /dev/null and b/examples/gcil/dataset/pubmed/0.01_1_6.npz differ
diff --git a/examples/gcil/dataset/pubmed/0.01_1_7.npz b/examples/gcil/dataset/pubmed/0.01_1_7.npz
new file mode 100644
index 00000000..6bf1a8b8
Binary files /dev/null and b/examples/gcil/dataset/pubmed/0.01_1_7.npz differ
diff --git a/examples/gcil/dataset/pubmed/0.01_1_8.npz b/examples/gcil/dataset/pubmed/0.01_1_8.npz
new file mode 100644
index 00000000..9f530d00
Binary files /dev/null and b/examples/gcil/dataset/pubmed/0.01_1_8.npz differ
diff --git a/examples/gcil/dataset/pubmed/0.01_1_9.npz b/examples/gcil/dataset/pubmed/0.01_1_9.npz
new file mode 100644
index 00000000..dd52acd0
Binary files /dev/null and b/examples/gcil/dataset/pubmed/0.01_1_9.npz differ
diff --git a/examples/gcil/gcil_trainer.py b/examples/gcil/gcil_trainer.py
new file mode 100644
index 00000000..54d5bca3
--- /dev/null
+++ b/examples/gcil/gcil_trainer.py
@@ -0,0 +1,425 @@
+import os
+
+# os.environ['CUDA_VISIBLE_DEVICES'] = '0'
+# os.environ['TL_BACKEND'] = 'torch'
+os.environ['TF_CPP_MIN_LOG_LEVEL'] = '2' 
+# 0:Output all; 1:Filter out INFO; 2:Filter out INFO and WARNING; 3:Filter out INFO, WARNING, and ERROR
+
+
+import numpy as np
+
+from gammagl.models import GCILModel, LogReg
+from gammagl.utils import  mask_to_index, add_self_loops
+from aug import random_aug
+from params import set_params
+
+from gammagl.data import Graph
+from gammagl.datasets import Planetoid
+
+from sklearn.metrics import f1_score
+import scipy.sparse as sp
+import pandas as pd
+import warnings
+import tensorlayerx as tlx
+from tensorlayerx.optimizers import Adam
+from tensorlayerx.model import TrainOneStep, WithLoss
+
+
+warnings.filterwarnings('ignore')
+
+class SemiSpvzLoss(WithLoss):
+    def __init__(self, net, loss_fn):
+        super(SemiSpvzLoss, self).__init__(backbone=net, loss_fn=loss_fn)
+
+    def forward(self, data, y): 
+        z1, z2, h1, h2 = self.backbone_network(data['graph1'], data['feat1'], data['attr1'], data['graph2'], data['feat2'], data['attr2'])
+        std_x = tlx.sqrt(h1.var(dim=0) + 0.0001)
+        std_y = tlx.sqrt(h2.var(dim=0) + 0.0001)
+
+        std_loss = tlx.ops.reduce_sum(tlx.sqrt((1 - std_x)**2)) / 2 + tlx.ops.reduce_sum(tlx.sqrt((1 - std_y)**2)) / 2
+        
+        c = tlx.matmul(z1.T, z2) / data['num_nodes']
+        c1 = tlx.matmul(z1.T, z1) / data['num_nodes']
+        c2 = tlx.matmul(z2.T, z2) / data['num_nodes']
+
+        loss_inv = -tlx.diag(c).sum() 
+        iden = tlx.convert_to_tensor(np.eye(c.shape[0])).to(args.device)
+        loss_dec1 = (iden - c1).pow(2).sum()
+        loss_dec2 = (iden - c2).pow(2).sum()
+        
+        loss = data['alpha']*loss_inv + data['beta'] * \
+            (loss_dec1 + loss_dec2) + data['gamma']*std_loss
+        return loss
+
+class LogRegLoss(WithLoss):
+    def __init__(self, net, loss_fn):
+        super(LogRegLoss, self).__init__(backbone=net, loss_fn=loss_fn)
+
+    def forward(self, data, y): 
+        logits = self.backbone_network(data['train_embs'])
+        loss = self._loss_fn(logits, y)
+        return loss
+
+def calculate_acc(logits, y, metrics):
+    metrics.update(logits, y) 
+    rst = metrics.result() 
+    metrics.reset() 
+    return rst
+
+def normalize_adj(adj):
+    """Symmetrically normalize adjacency matrix."""
+    adj = sp.coo_matrix(adj)
+    rowsum = np.array(np.abs(adj.A).sum(1))
+    d_inv_sqrt = np.power(rowsum, -0.5).flatten()
+    d_inv_sqrt[np.isinf(d_inv_sqrt)] = 0.
+    d_mat_inv_sqrt = sp.diags(d_inv_sqrt)
+    return adj.dot(d_mat_inv_sqrt).transpose().dot(d_mat_inv_sqrt).tocoo()
+
+def get_dataset(name):
+    if name == 'cora':
+        dataset = Planetoid(root="", name='cora')
+    elif name == 'pubmed':
+        dataset = Planetoid(root="", name='pubmed')
+    else:
+        raise ValueError(f"Unknown dataset name: {name}")
+    graph = dataset[0]
+
+    num_nodes = graph.num_nodes
+    edge_index = graph.edge_index
+    scope = add_self_loops(edge_index)[0]
+    adj = sp.coo_matrix((np.ones(tlx.get_tensor_shape(edge_index)[1]),
+                          (tlx.to_device(edge_index[0, :], "cpu"), tlx.to_device(edge_index[1, :], "cpu"))),
+                          shape=[num_nodes, num_nodes])
+
+
+    feat = graph.x
+    feat = tlx.convert_to_tensor(feat, dtype=tlx.float32)
+
+
+    label = graph.y
+    num_class = label.max()+1
+
+    idx_train = mask_to_index(graph.train_mask)
+    idx_test = mask_to_index(graph.test_mask)
+    idx_val = mask_to_index(graph.val_mask)
+
+    laplace = sp.eye(adj.shape[0]) - normalize_adj(adj)
+    
+    return adj, feat, label, num_class, idx_train, idx_val, idx_test, laplace, scope
+
+
+def test(embeds, labels, num_class, train_idx, val_idx, test_idx):
+
+    train_embs = embeds[train_idx]
+    val_embs = embeds[val_idx]
+    test_embs = embeds[test_idx]
+
+    label = labels.to(args.device)
+
+    train_labels = label[train_idx]
+    val_labels = label[val_idx]
+    test_labels = label[test_idx]
+
+    ''' Linear Evaluation '''
+    # print(train_embs.shape)
+    logreg = LogReg(train_embs.shape[1], num_class)
+    optimizer = Adam(lr=args.lr2, weight_decay=args.wd2)
+
+    # logreg = logreg.to(args.device)
+    train_weights = logreg.trainable_weights
+    loss_func = LogRegLoss(logreg, tlx.losses.softmax_cross_entropy_with_logits)
+    
+    train_one_step = TrainOneStep(loss_func, optimizer, train_weights)
+
+    best_val_acc = 0
+    eval_acc = 0
+    data = {
+        'train_embs': train_embs
+    }
+    for epoch in range(800):
+        logreg.train()
+        train_one_step(data, train_labels)
+        
+        logreg.eval()
+        logits = logreg(train_embs)
+        preds = tlx.argmax(logits, axis=1)
+        train_acc = tlx.reduce_sum(preds == train_labels).float() / train_labels.shape[0]
+
+        val_logits = logreg(val_embs)
+        test_logits = logreg(test_embs)
+
+        val_preds = tlx.argmax(val_logits, axis=1)
+        test_preds = tlx.argmax(test_logits, axis=1)
+
+        val_acc = tlx.reduce_sum(val_preds == val_labels).float() / val_labels.shape[0]
+        test_acc = tlx.reduce_sum(test_preds == test_labels).float() / test_labels.shape[0]
+
+        test_f1_macro = f1_score(
+            test_labels.cpu(), test_preds.cpu(), average='macro')
+        test_f1_micro = f1_score(
+            test_labels.cpu(), test_preds.cpu(), average='micro')
+        if val_acc >= best_val_acc:
+            best_val_acc = val_acc
+            if test_acc > eval_acc:
+                test_f1_macro_ll = test_f1_macro
+                test_f1_micro_ll = test_f1_micro
+
+    print('Epoch:{}, train_acc:{:.4f}, Macro:{:4f}, Micro:{:4f}'.format(epoch, train_acc, test_f1_macro_ll,
+                                                                        test_f1_micro_ll))
+
+    return test_f1_macro_ll, test_f1_micro_ll
+
+
+def train(params):
+    path = "./dataset/" + args.dataset
+    adj, feat, labels, num_class, train_idx, val_idx, test_idx, laplace, scope = get_dataset(args.dataset)
+
+    adj = adj + sp.eye(adj.shape[0])
+
+    edge_index = np.vstack((adj.nonzero()[0], adj.nonzero()[1]))
+    x = tlx.convert_to_tensor(feat, dtype=tlx.float32)  
+    graph = Graph(x=x, edge_index=tlx.convert_to_tensor(edge_index, dtype=tlx.int64))
+    
+    if args.dataset == 'pubmed':
+        new_adjs = []
+        for i in range(10):
+            new_adjs.append(sp.load_npz(path + "/0.01_1_" + str(i) + ".npz"))
+        adj_num = len(new_adjs)
+        adj_inter = int(adj_num / args.num)
+        sele_adjs = []
+        for i in range(args.num + 1):
+            try:
+                if i == 0:
+                    sele_adjs.append(new_adjs[i])
+                else:
+                    sele_adjs.append(new_adjs[i * adj_inter - 1])
+            except IndexError:
+                pass
+        print("Number of select adjs:", len(sele_adjs))
+        epoch_inter = args.epoch_inter
+    elif args.dataset == 'cora':
+        sele_adjs = []
+        for i in range(7):
+            sele_adjs.append(sp.load_npz(path + "/0.01_1_" + str(i) + ".npz"))
+        epoch_inter = args.epoch_inter
+
+    if args.gpu != -1:
+        args.device = 'cuda:{}'.format(args.gpu) 
+    else:
+        args.device = 'cpu'
+
+    in_dim = feat.shape[1]
+    
+    model = GCILModel(in_dim, args.hid_dim, args.out_dim, args.n_layers, args.use_mlp)
+    # model = model.to(args.device)
+
+    optimizer = tlx.optimizers.Adam(
+    lr=args.lr1,
+    weight_decay=args.wd1
+    )
+    train_weights = model.trainable_weights
+    loss_func = SemiSpvzLoss(model, tlx.losses.softmax_cross_entropy_with_logits)
+    
+    train_one_step = TrainOneStep(loss_func, optimizer, train_weights)
+
+    #### SpCo ######
+    theta = 1
+    delta = np.ones(adj.shape) * args.delta_origin
+    delta_add = delta
+    delta_dele = delta
+    num_node = adj.shape[0]
+    range_node = np.arange(num_node)
+    ori_graph = graph
+    new_graph = ori_graph
+
+    new_adj = adj.tocsc()
+    ori_attr = tlx.convert_to_tensor(new_adj[new_adj.nonzero()])[0]  
+    ori_diag_attr = tlx.convert_to_tensor(new_adj[range_node, range_node])[0] 
+    new_attr = tlx.convert_to_tensor(new_adj[new_adj.nonzero()])[0]
+    new_diag_attr = tlx.convert_to_tensor(new_adj[range_node, range_node])[0]
+    j = 0
+    
+    # Training Loop
+    for epoch in range(params['epoch']):
+        model.train()
+        # train_loss = train_one_step(data, graph.y)
+
+        graph1_, attr1, feat1 = random_aug(new_graph, new_attr, new_diag_attr, feat, args.dfr, args.der, args.device)
+        graph2_, attr2, feat2 = random_aug(ori_graph, ori_attr, ori_diag_attr, feat, args.dfr, args.der, args.device)
+
+        graph1 = Graph(x=graph1_.x.to(args.device), edge_index=graph1_.edge_index.to(args.device))
+        graph2 = Graph(x=graph2_.x.to(args.device), edge_index=graph2_.edge_index.to(args.device))
+
+        attr1 = attr1.to(args.device)
+        attr2 = attr2.to(args.device)
+
+        feat1 = feat1.to(args.device)
+        feat2 = feat2.to(args.device)
+
+        data = {
+            "x": graph.x,
+            "y": graph.y,
+            "edge_index": edge_index,
+            "in_dim": in_dim,
+            "num_nodes": graph.num_nodes,
+            "graph1": graph1,
+            "graph2": graph2,
+            "attr1": attr1,
+            "attr2": attr2,
+            "feat1": feat1,
+            "feat2": feat2,
+            "alpha": params['alpha'],
+            'beta': params['beta'],
+            "gamma": params['gamma']
+        }
+        
+        loss = train_one_step(data, graph.y)
+
+        z1, z2, h1, h2 = model(graph1, feat1, attr1, graph2, feat2, attr2)
+        c = tlx.matmul(z1.T, z2) / graph.num_nodes
+        c1 = tlx.matmul(z1.T, z1) / graph.num_nodes
+        c2 = tlx.matmul(z2.T, z2) / graph.num_nodes
+        iden = tlx.convert_to_tensor(np.eye(c.shape[0])).to(args.device)
+        # loss.backward()
+
+        # Print the results
+        print('Epoch={:03d}, loss={:.4f}, loss_inv={:.4f}, loss_dec={:.4f}'.format(
+        epoch, 
+        loss.item(),  
+        -tlx.diag(c).sum().item(),  
+        (iden - c1).pow(2).sum().item() + (iden - c2).pow(2).sum().item()  # 同样处理loss_dec
+        ))
+
+
+        if args.dataset == 'pubmed':
+            if (epoch - 1) % epoch_inter == 0:
+                try:
+                    print("================================================")
+                    delta = args.lam * sele_adjs[int(epoch / epoch_inter)]
+                    new_adj = adj + delta
+
+                    nonzero_indices = np.array(new_adj.nonzero())
+                    edge_index = tlx.ops.convert_to_tensor(nonzero_indices, dtype=tlx.int64)
+                    
+                    x = tlx.convert_to_tensor(feat, dtype=tlx.float32)
+                    new_graph = Graph(x=x, edge_index=edge_index)
+                    
+                    new_attr = tlx.convert_to_tensor(new_adj[new_adj.nonzero()], dtype=tlx.float32)[0]  # 使用 tlx.float32
+                    new_diag_attr = tlx.convert_to_tensor(new_adj[range_node, range_node], dtype=tlx.float32)[0]  # 使用 tlx.float32
+
+                except IndexError:
+                    pass
+
+
+        elif args.dataset == 'cora':
+            flag = (epoch - 1) % epoch_inter
+            if flag == 0:
+                try:
+                    print("================================================")
+                    delta = args.lam * sele_adjs[(epoch - 1)//epoch_inter]
+                    new_adj = adj + delta
+                    
+                    nonzero_indices = np.array(new_adj.nonzero())
+                    edge_index = tlx.ops.convert_to_tensor(nonzero_indices, dtype=tlx.int64)
+
+                    x = tlx.convert_to_tensor(feat, dtype=tlx.float32)
+                    new_graph = Graph(x=x, edge_index=edge_index)
+
+                    new_attr = tlx.convert_to_tensor(new_adj[new_adj.nonzero()], dtype=tlx.float32)[0]  # 使用 tlx.float32
+                    new_diag_attr = tlx.convert_to_tensor(new_adj[range_node, range_node], dtype=tlx.float32)[0]  # 使用 tlx.float32
+                except IndexError:
+                    pass
+
+
+    # Final Testing
+    # Remove direct use of graph.to(args.device) and replace it with moving the graph attributes manually
+    graph.x = graph.x.to(args.device)  # Move features to device
+    graph.edge_index = graph.edge_index.to(args.device)  # Move edge_index to device
+    # Manually add self-loops without using remove_self_loop and add_self_loop methods
+    num_nodes = graph.num_nodes
+    edge_index = graph.edge_index.cpu().numpy()  # Convert to numpy for edge_index
+    self_loops = np.vstack([np.arange(num_nodes), np.arange(num_nodes)])
+    edge_index = np.hstack([edge_index, self_loops])  # Add self loops to edge_index
+    graph.edge_index = tlx.convert_to_tensor(edge_index, dtype=tlx.int64, device=args.device)  # Convert back to TensorLayerX tensor
+
+    feat = feat.to(args.device)
+
+    # Manually reconstruct adjacency matrix from edge_index
+    new_adj = sp.coo_matrix((np.ones(graph.edge_index.shape[1]), 
+                            (graph.edge_index[0].cpu().numpy(), graph.edge_index[1].cpu().numpy())),
+                            shape=(num_nodes, num_nodes)).tocsc()
+
+    # Reconstruct attributes from adjacency matrix
+    attr = tlx.convert_to_tensor(new_adj[new_adj.nonzero()], dtype=tlx.float32).to(args.device)
+
+    # Get embeddings using the model
+    embeds = model.get_embedding(graph, feat, attr)
+
+    # Evaluation
+    test_f1_macro_ll = 0
+    test_f1_micro_ll = 0
+    macros = []
+    micros = []
+
+    # Run the test function
+    test_f1_macro_ll, test_f1_micro_ll = test(embeds, labels, num_class, train_idx, val_idx, test_idx)
+
+    # Append results to macros and micros
+    macros.append(test_f1_macro_ll)
+    micros.append(test_f1_micro_ll)
+
+    # Convert lists to tensors
+    macros = tlx.convert_to_tensor(macros, dtype=tlx.float32)
+    micros = tlx.convert_to_tensor(micros, dtype=tlx.float32)
+
+    # Convert arguments to dictionary (if necessary)
+    config = vars(args)
+
+
+    config['test_f1_macro'] = test_f1_macro_ll
+    config['test_f1_micro'] = test_f1_micro_ll
+    print(config)
+
+    return tlx.ops.reduce_mean(macros).cpu().numpy().item(), tlx.ops.reduce_mean(micros).cpu().numpy().item()
+
+
+
+def main(args):
+    print(args)
+    macros, micros = [], []
+
+    df = pd.DataFrame()
+    params = {}
+    params['alpha'], params['beta'], params['gamma'], params['epoch'] = 1, 0.015, 0, 100
+
+    for i in range(10):
+        ma, mi = train(params)
+        macros.append(ma)
+        micros.append(mi)
+        print(ma, mi)
+        print(params)
+
+    print(params)
+
+    micros = tlx.convert_to_tensor(micros)
+    macros = tlx.convert_to_tensor(macros)
+
+    print('AVG Micro:{:.4f}, Std:{:.4f}, Macro:{:.4f}, Std:{:.4f}'.format(
+    tlx.ops.reduce_mean(micros).cpu().numpy().item(),
+    tlx.ops.reduce_std(micros).cpu().numpy().item(),
+    tlx.ops.reduce_mean(macros).cpu().numpy().item(),
+    tlx.ops.reduce_std(macros).cpu().numpy().item())
+)
+
+
+
+
+if __name__ == '__main__':
+    args = set_params()
+
+    if args.gpu >= 0:
+        tlx.set_device("GPU", args.gpu)
+    else:
+        tlx.set_device("CPU")
+
+    main(args)
diff --git a/examples/gcil/params.py b/examples/gcil/params.py
new file mode 100644
index 00000000..e60ce024
--- /dev/null
+++ b/examples/gcil/params.py
@@ -0,0 +1,102 @@
+import argparse
+import sys
+
+argv = sys.argv
+dataset = argv[1]
+
+if len(argv) < 2:
+    print("错误: 缺少必需的 '--dataset' 参数。")
+    sys.exit(1)
+
+def cora_params():
+    parser = argparse.ArgumentParser()
+
+    parser.add_argument('--dataset', type=str, default="cora")    
+    parser.add_argument('--gpu', type=int, default=0)
+    parser.add_argument('--lr1', type=float, default=1e-3)
+    parser.add_argument('--lr2', type=float, default=1e-2)
+    parser.add_argument('--wd1', type=float, default=0)
+    parser.add_argument('--wd2', type=float, default=1e-4)
+    parser.add_argument('--lambd', type=float, default=1e-3)
+    parser.add_argument('--der', type=float, default=0.4)
+    parser.add_argument('--dfr', type=float, default=0.1) 
+    parser.add_argument('--use_mlp', action='store_true', default=False)
+    parser.add_argument('--hid_dim', type=int, default=512)  
+    parser.add_argument('--out_dim', type=int, default=512) 
+    parser.add_argument('--n_layers', type=int, default=2)
+
+
+    #####################################
+    ## new parameters
+    parser.add_argument('--delta_origin', type=float, default=0.5)
+    parser.add_argument('--sin_iter', type=int, default=3) 
+
+    parser.add_argument('--lam', type=float, default=0.5) 
+    parser.add_argument('--epochs', type=int, default=40)
+    parser.add_argument('--turn', type=int, default=15)
+    parser.add_argument('--epsilon', type=float, default=0.01)
+    parser.add_argument('--scope_flag', type=int, default=1)
+    parser.add_argument('--epoch_inter', type=int, default=15)
+    #####################################
+
+    args, _ = parser.parse_known_args()
+    return args
+
+def pubmed_params():
+    parser = argparse.ArgumentParser()
+
+    parser.add_argument('--dataset', type=str, default="pubmed")    
+    parser.add_argument('--gpu', type=int, default=0) 
+    parser.add_argument('--lr1', type=float, default=1e-3)
+    parser.add_argument('--lr2', type=float, default=1e-2)
+    parser.add_argument('--wd1', type=float, default=0)
+    parser.add_argument('--wd2', type=float, default=1e-4)
+    parser.add_argument('--lambd', type=float, default=1e-3)
+    parser.add_argument('--der', type=float, default=0.5)
+    parser.add_argument('--dfr', type=float, default=0.3) 
+    parser.add_argument('--use_mlp', action='store_true', default=False)
+    parser.add_argument('--hid_dim', type=int, default=512)  
+    parser.add_argument('--out_dim', type=int, default=512) 
+    parser.add_argument('--n_layers', type=int, default=2)
+    parser.add_argument('--seed', type=int, default=2)
+
+
+    #####################################
+    ## new parameters
+    parser.add_argument('--delta_origin', type=float, default=0.5)
+    parser.add_argument('--theta', type=float, default=1.)  
+    parser.add_argument('--sin_iter', type=int, default=1)
+
+    # 20 / 10 labeled nodes per class
+    parser.add_argument('--lam', type=float, default=0.1)
+    parser.add_argument('--epochs', type=int, default=170)
+    parser.add_argument('--num', type=int, default=7)
+    parser.add_argument('--epoch_inter', type=int, default=15)
+    parser.add_argument('--scope_flag', type=int, default=1)
+
+    '''
+    # 5 labeled nodes per class
+    parser.add_argument('--lam', type=float, default=0.5)
+    parser.add_argument('--epochs', type=int, default=150)
+    parser.add_argument('--num', type=int, default=2)
+    parser.add_argument('--epoch_inter', type=int, default=15)
+    parser.add_argument('--scope_flag', type=int, default=1)
+    '''
+    #####################################
+
+    args, _ = parser.parse_known_args()
+    return args
+
+def set_params():
+    if dataset == "cora":
+        args = cora_params()
+    elif dataset == "pubmed":
+        args = pubmed_params()
+    else:
+        # 处理无效 dataset 的情况
+        print(f"错误: 不支持的数据集 '{dataset}'")
+        sys.exit(1)  # 退出程序，返回错误代码
+    return args
+
+# 主程序中调用
+args = set_params()
diff --git a/gammagl/models/__init__.py b/gammagl/models/__init__.py
index 24ed6176..94862783 100644
--- a/gammagl/models/__init__.py
+++ b/gammagl/models/__init__.py
@@ -62,6 +62,7 @@
 from .dna import DNAModel
 from .dfad import DFADModel, DFADGenerator
 from .fatragnn import FatraGNNModel, Graph_Editer
+from .gcil import GCILModel, LogReg 
 
 __all__ = [
     'HeCo',
@@ -130,7 +131,9 @@
     'DNAModel',
     'DFADModel',
     'DFADGenerator',
-    'FatraGNNModel'
+    'FatraGNNModel',
+    'GCILModel',
+    'LogReg'
 ]
 
 classes = __all__
diff --git a/gammagl/models/gcil.py b/gammagl/models/gcil.py
new file mode 100644
index 00000000..13d08422
--- /dev/null
+++ b/gammagl/models/gcil.py
@@ -0,0 +1,91 @@
+import tensorlayerx.nn as nn
+import tensorlayerx as tlx
+from gammagl.models.gcn import GCNModel as GCN
+
+
+class GCILModel(nn.Module):
+    r"""GCIL Model proposed in '"Graph Contrastive Invariant Learning from the Causal Perspective"
+        <https://arxiv.org/pdf/2401.12564v2.pdf>'_ paper.
+            
+        Parameters
+        ----------
+        in_dim: int
+            Input feature dimension.
+        hid_dim: int
+            Hidden layer dimension.
+        out_dim: int
+            Output feature dimension, usually the number of classes or the desired embedding size.
+        n_layers: int
+            Number of layers for the GCN model (if GCN is used).
+        use_mlp: bool, optional
+            If True, the model will use an MLP backbone instead of GCN. Default is False.
+    """
+    def __init__(self, in_dim, hid_dim, out_dim, n_layers, use_mlp=False, drop_rate=0.2):
+        super().__init__()
+        if not use_mlp:
+            self.backbone = GCN(in_dim, hid_dim, out_dim, drop_rate=drop_rate, num_layers=n_layers)
+        else:
+            self.backbone = MLP(in_dim, hid_dim, out_dim)
+            
+
+    def get_embedding(self, graph, feat, attr):
+        x = graph.x
+        num_nodes = graph.num_nodes
+        edge_index = graph.edge_index
+        edge_weight = graph.edge_weight if hasattr(graph, 'edge_weight') else None
+        out = self.backbone(x, edge_index, edge_weight, num_nodes)
+        return out.detach()
+
+    def forward(self, graph1, feat1, attr1, graph2, feat2, attr2):
+        x1 = graph1.x  
+        x2 = graph2.x 
+        
+        num_nodes1 = graph1.num_nodes
+        num_nodes2 = graph2.num_nodes
+        
+        edge_index1 = graph1.edge_index 
+        edge_index2 = graph2.edge_index 
+        
+        edge_weight1 = graph1.edge_weight if hasattr(graph1, 'edge_weight') else None
+        edge_weight2 = graph2.edge_weight if hasattr(graph2, 'edge_weight') else None
+        
+
+        h1 = self.backbone(x1, edge_index1, edge_weight1, num_nodes1)
+        h2 = self.backbone(x2, edge_index2, edge_weight2, num_nodes2)
+
+        z1 = (h1 - h1.mean(0)) / h1.std(0)
+        z2 = (h2 - h2.mean(0)) / h2.std(0)
+
+        return z1, z2, h1, h2
+
+class LogReg(nn.Module):
+    def __init__(self, hid_dim, out_dim):
+        super(LogReg, self).__init__()
+        # self.fc = nn.Linear(hid_dim, out_dim)
+        self.fc = nn.Linear(in_features=hid_dim, out_features=out_dim)
+
+    def forward(self, x):
+        z = self.fc(x)
+        return z
+
+
+class MLP(nn.Module):
+    def __init__(self, nfeat, nhid, nclass, use_bn=True):
+        super(MLP, self).__init__()
+
+        self.layer1 = nn.Linear(in_features=nfeat, out_features=nhid)
+        self.layer2 = nn.Linear(in_features=nhid, out_features=nclass)
+
+        self.bn = nn.BatchNorm1d(nhid)
+        self.use_bn = use_bn
+        self.act_fn = nn.ReLU()
+
+    def forward(self, _, x):
+        x = self.layer1(x)
+        if self.use_bn:
+            x = self.bn(x)
+
+        x = self.act_fn(x)
+        x = self.layer2(x)
+
+        return x
\ No newline at end of file