{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.7.10","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"gpu","dataSources":[{"sourceId":10418,"databundleVersionId":862236,"sourceType":"competition"},{"sourceId":252838,"sourceType":"datasetVersion","datasetId":106096}],"dockerImageVersionId":30097,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Researchera - 计算机视觉项目\n# Human Protein Atlas Image Classification\n\n在此Kaggle竞赛中，参赛者需要建立机器学习模型实现显微镜下的混合蛋白质图像的分类（即一张图片拥有多个标签，为多标签分类问题）。[The Human Protein Atlas](https://www.proteinatlas.org/) 会使用这些模型来创建智能显微系统，进而从高分辨率图像中准确识别蛋白质类型。\n\n蛋白质是人体细胞功能的执行者，维持着生命的正常运转。在之前的研究中，蛋白质图像分类的模式较为单一，图片上通常只有一种或少部分类别细胞。但是为了能够深入理解人类细胞机制的复杂性，此竞赛要求参赛者在存在多个不同类别细胞的图片中准确识别蛋白质类型。\n\n细胞蛋白质的可视化常用于生物学研究，并且可以推动医药研究的突破。目前，由于高清显微镜技术的快速发展，图像数据产生的速度已经超过了数据处理的速度。因此，主办方希望通过此次竞赛研究生物图像自动分析的技术，以此推动人类对细胞和疾病的理解。\n\n竞赛优胜者将受邀与 The Human Protein Atlas 团队一起在 Nature 旗下子刊 [*Nature Methods*](https://www.nature.com/nmeth) 上发文，此期刊最新影响因子高达28.467！\n\n>**提示：**Code 和 Markdown 区域可通过 **Shift + Enter** 快捷键运行。此外，Markdown可以通过双击进入编辑模式。\n\n我们将这个notebook分为不同的步骤，你可以使用下面的链接来浏览此notebook。\n\n* [Step 1](#step1): 导入数据\n* [Step 2](#step2): 数据研究\n* [Step 3](#step3): 构建PyTorch数据集\n* [Step 4](#step4): 搭建四通道CNN模型\n* [Step 5](#step5): Focal Loss\n* [Step 6](#step6): 模型训练\n* [Step 7](#step7): 后处理与提交\n\n在该项目中包含了如下的问题：\n\n* [问题 1](#question1): 查阅资料并回答：如何解决训练数据中标签分布不均衡的问题？\n* [问题 2](#question2): 此处可否简单采用`sklearn.model_selection`中的`StratifiedKFold`方法进行数据集切分？\n* [问题 3](#question3): 简要概述Focal Loss与BCE Loss的差别。","metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5"}},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport random\nimport os\nimport cv2\nimport gc\nimport warnings\nfrom sklearn.metrics import f1_score\nfrom sklearn.exceptions import UndefinedMetricWarning\nimport scipy.optimize as opt\nfrom collections import defaultdict, Counter\nfrom tqdm.notebook import tqdm\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport torch\nfrom torch import nn, optim\nfrom torch.nn import functional as F\nfrom torch.utils.data import Dataset, DataLoader\nfrom torchvision.models import resnet50\n%matplotlib inline","metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","execution":{"iopub.status.busy":"2024-10-14T11:32:06.218633Z","iopub.execute_input":"2024-10-14T11:32:06.219022Z","iopub.status.idle":"2024-10-14T11:32:06.231831Z","shell.execute_reply.started":"2024-10-14T11:32:06.218989Z","shell.execute_reply":"2024-10-14T11:32:06.230692Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<a id='step1'></a>\n# 1. 导入数据\n\n在此竞赛中，主办方给出了124288个训练样本和46808个测试样本。每个样本均由四张512$\\times$512像素的单通道图片组成，分别代表红色（Red）、绿色（Green）、蓝色（Blue）和黄色（Yellow）四个通道。图片文件名由样本代号+通道代号组成，例如，样本代号为`00070df0-bbc3-11e8-b2bc-ac1f6b6435d0`的样本一共对应四张图片，分别为：\n* 00070df0-bbc3-11e8-b2bc-ac1f6b6435d0_red.png\n* 00070df0-bbc3-11e8-b2bc-ac1f6b6435d0_green.png\n* 00070df0-bbc3-11e8-b2bc-ac1f6b6435d0_blue.png\n* 00070df0-bbc3-11e8-b2bc-ac1f6b6435d0_yellow.png\n\n训练集中所有样本的代号记录在`train.csv`中，测试集中所有样本的代号记录在`sample_submission.csv`中。数据集中共含有28种人类蛋白质，其中训练集的标签在`train.csv`中给出。下面我们分别导入：","metadata":{}},{"cell_type":"code","source":"# 文件路径\nTRAIN = '../input/human-protein-atlas-image-classification/train/'\nTEST = '../input/human-protein-atlas-image-classification/test/'\nLABELS = '../input/human-protein-atlas-image-classification/train.csv'\nSUBMIT = '../input/human-protein-atlas-image-classification/sample_submission.csv'","metadata":{"execution":{"iopub.status.busy":"2024-10-14T11:32:06.233290Z","iopub.execute_input":"2024-10-14T11:32:06.233616Z","iopub.status.idle":"2024-10-14T11:32:06.245486Z","shell.execute_reply.started":"2024-10-14T11:32:06.233586Z","shell.execute_reply":"2024-10-14T11:32:06.244570Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# 标签中的数字与蛋白质种类的对应关系\nname_label_dict = {\n0:  'Nucleoplasm',\n1:  'Nuclear membrane',\n2:  'Nucleoli',   \n3:  'Nucleoli fibrillar center',\n4:  'Nuclear speckles',\n5:  'Nuclear bodies',\n6:  'Endoplasmic reticulum',   \n7:  'Golgi apparatus',\n8:  'Peroxisomes',\n9:  'Endosomes',\n10:  'Lysosomes',\n11:  'Intermediate filaments',\n12:  'Actin filaments',\n13:  'Focal adhesion sites',   \n14:  'Microtubules',\n15:  'Microtubule ends',  \n16:  'Cytokinetic bridge',   \n17:  'Mitotic spindle',\n18:  'Microtubule organizing center',  \n19:  'Centrosome',\n20:  'Lipid droplets',\n21:  'Plasma membrane',   \n22:  'Cell junctions', \n23:  'Mitochondria',\n24:  'Aggresome',\n25:  'Cytosol',\n26:  'Cytoplasmic bodies',   \n27:  'Rods & rings' }\n\nN_CLASSES = len(name_label_dict)\nprint('The number of classes of Human Proteins in this competition is {}'.format(N_CLASSES))","metadata":{"execution":{"iopub.status.busy":"2024-10-14T11:32:06.247220Z","iopub.execute_input":"2024-10-14T11:32:06.247577Z","iopub.status.idle":"2024-10-14T11:32:06.258134Z","shell.execute_reply.started":"2024-10-14T11:32:06.247531Z","shell.execute_reply":"2024-10-14T11:32:06.257235Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# 训练集\ndf = pd.read_csv(LABELS)\ndf.head()","metadata":{"execution":{"iopub.status.busy":"2024-10-14T11:32:06.259479Z","iopub.execute_input":"2024-10-14T11:32:06.259779Z","iopub.status.idle":"2024-10-14T11:32:06.311017Z","shell.execute_reply.started":"2024-10-14T11:32:06.259742Z","shell.execute_reply":"2024-10-14T11:32:06.310111Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"可以看出一张图片上可能存在多个种类的蛋白质，即这个竞赛的任务实际上是一个多标签分类（Multi-Label Classification）问题。","metadata":{}},{"cell_type":"code","source":"# 测试集\nsub_df = pd.read_csv(SUBMIT)\nsub_df.head()","metadata":{"execution":{"iopub.status.busy":"2024-10-14T11:32:06.312341Z","iopub.execute_input":"2024-10-14T11:32:06.312697Z","iopub.status.idle":"2024-10-14T11:32:06.335471Z","shell.execute_reply.started":"2024-10-14T11:32:06.312659Z","shell.execute_reply":"2024-10-14T11:32:06.334707Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<a id='step2'></a>\n# 2. 数据研究\n\n在划分训练集与验证集时，要保证两个数据集中的标签分布尽可能一致，这样才能通过验证集有效评估模型的泛化能力。因此，在划分数据集之前，需要首先研究数据标签的分布情况。","metadata":{}},{"cell_type":"code","source":"cls_counts = Counter(cls for classes in df['Target'].str.split() for cls in classes)\ncounts_x = [i[1] for i in cls_counts.most_common(N_CLASSES)]\ncounts_y = [name_label_dict[int(i[0])] for i in cls_counts.most_common(N_CLASSES)]\nplt.figure(figsize=(8,8))\nsns.barplot(y=counts_y, x=counts_x)","metadata":{"execution":{"iopub.status.busy":"2024-10-14T11:32:06.336465Z","iopub.execute_input":"2024-10-14T11:32:06.336733Z","iopub.status.idle":"2024-10-14T11:32:06.811519Z","shell.execute_reply.started":"2024-10-14T11:32:06.336708Z","shell.execute_reply":"2024-10-14T11:32:06.810669Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"不难看出，`Nucleoplasm`和`Cytosl`两种蛋白质占据了标签的绝大部分，而`Lysosomes`、`Microtubule ends`和`Rods & rings`这三类蛋白质数量却相当少。经统计可得：数量最多的`Nucleoplasm`占比高达 25%，而数量最少的`Rods & rings`占比只有 0.02%。因此，此竞赛所处理的数据是一个典型的非均衡数据集。\n\n<a id='question1'></a>\n### __问题 1:__\n\n查阅资料并回答：如何解决训练数据中标签分布不均衡的问题？\n\n__回答:__ ","metadata":{}},{"cell_type":"markdown","source":"1.对占比较多的样本进行欠采样，对占比较少的样本进行过采样。2.在模型参数中加入class_weight并增加占比较少的样本的权重。","metadata":{}},{"cell_type":"markdown","source":"## 2.1 StratifiedKFold\n\n为保证训练集与验证集中的标签分布的一致性，我们使用`StratifiedKFold`方法进行数据集的划分。`StratifiedKFold`为分层采样交叉切分，确保训练集、测试集中各类别样本的比例与原始数据集中相同。","metadata":{}},{"cell_type":"code","source":"# Fold划分\nn_folds = 10\nfold_cls_counts = defaultdict(int)\nfolds = [-1] * len(df)\nfor item in tqdm(df.sample(frac=1, random_state=42).itertuples(),total=len(df)):\n    cls = min(item.Target.split(), key=lambda cls: cls_counts[cls])\n    fold_counts = [(f, fold_cls_counts[f, cls]) for f in range(n_folds)]\n    min_count = min([count for _, count in fold_counts])\n    random.seed(item.Index)\n    fold = random.choice([f for f, count in fold_counts if count == min_count])\n    folds[item.Index] = fold\n    for cls in item.Target.split():\n        fold_cls_counts[fold, cls] += 1\ndf['fold'] = folds","metadata":{"execution":{"iopub.status.busy":"2024-10-14T11:32:06.812672Z","iopub.execute_input":"2024-10-14T11:32:06.812944Z","iopub.status.idle":"2024-10-14T11:32:07.623404Z","shell.execute_reply.started":"2024-10-14T11:32:06.812918Z","shell.execute_reply":"2024-10-14T11:32:07.622439Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df.head(10)","metadata":{"execution":{"iopub.status.busy":"2024-10-14T11:32:07.624633Z","iopub.execute_input":"2024-10-14T11:32:07.624962Z","iopub.status.idle":"2024-10-14T11:32:07.636028Z","shell.execute_reply.started":"2024-10-14T11:32:07.624932Z","shell.execute_reply":"2024-10-14T11:32:07.635151Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"可以看到，每个样本都被划入相应的Fold中。接下来，我们只需要指定用作验证集的Fold，便可以完成训练集与验证集的划分。\n\n<a id='question2'></a>\n### __问题 2:__\n\n此处可否简单采用`sklearn.model_selection`中的`StratifiedKFold`方法进行数据集切分？\n\n__回答:__ \n","metadata":{}},{"cell_type":"markdown","source":"不可以。因为可以看到我们的样本是一个多标签分类的样本，sklearn中的StratifiedKFold工具一般是用于单标签分类的，它无法划分使得同时满足所有类在每个folder中的分布与原样本相同，所以需要我们手动实现。在手动实现时，我们规定每个样本中出现频率最小的标签作为分类标签，以此满足标签分布均衡化的要求。","metadata":{}},{"cell_type":"code","source":"valid_idx = 0\ntrain_df = df[df['fold']!=valid_idx][['Id', 'Target']].reset_index(drop=True)\nvalid_df = df[df['fold']==valid_idx][['Id', 'Target']].reset_index(drop=True)\nmulti = [1,1,1,1,1,1,1,1,4,4,4,1,1,1,1,4,\n        1,1,1,1,2,1,1,1,1,1,1,4]\n# 重复少见类样本\noversampled_data = []\nfor idx, row in train_df.iterrows():\n    targets = list(map(int, row['Target'].split()))\n    max_multi = max(multi[target] for target in targets)\n    for _ in range(max_multi):\n        oversampled_data.append(row)\n\ntrain_df = pd.DataFrame(oversampled_data).reset_index(drop=True)\nprint('There are {} samples in the training set.'.format(len(train_df)))\nprint('There are {} samples in the validation set.'.format(len(valid_df)))","metadata":{"execution":{"iopub.status.busy":"2024-10-14T11:32:07.638881Z","iopub.execute_input":"2024-10-14T11:32:07.639163Z","iopub.status.idle":"2024-10-14T11:32:12.118739Z","shell.execute_reply.started":"2024-10-14T11:32:07.639137Z","shell.execute_reply":"2024-10-14T11:32:12.117651Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 2.2 图片数据预览","metadata":{}},{"cell_type":"code","source":"# 定义RGBY图片读取函数\ndef open_rgby(path, id, size=None): #a function that reads RGBY image\n    colors = ['red','green','blue','yellow']\n    flags = cv2.IMREAD_GRAYSCALE\n    if size is None:\n        img = [cv2.imread(os.path.join(path, id+'_'+color+'.png'), flags).astype(np.float32)/255 \n               for color in colors]\n    else:\n        img = []\n        for color in colors:\n            src_img = cv2.imread(os.path.join(path, id+'_'+color+'.png'), flags)\n            tar_img = cv2.resize(src_img, (2*size, size), interpolation=cv2.INTER_CUBIC).astype(np.float32)/255\n            img.append(tar_img)\n\n    return np.stack(img, axis=0)","metadata":{"execution":{"iopub.status.busy":"2024-10-14T11:32:12.120479Z","iopub.execute_input":"2024-10-14T11:32:12.120875Z","iopub.status.idle":"2024-10-14T11:32:12.129340Z","shell.execute_reply.started":"2024-10-14T11:32:12.120841Z","shell.execute_reply":"2024-10-14T11:32:12.128452Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def make_image_row(image, subax, title):\n    subax[0].imshow(np.stack(image, axis=-1)[:,:,:3])\n    subax[0].set_title('Id: '+title)\n    subax[1].imshow(image[0], cmap=\"Reds\")\n    subax[1].set_title(\"Red Channel\")\n    subax[2].imshow(image[1], cmap=\"Greens\")\n    subax[2].set_title(\"Green Channel\")\n    subax[3].imshow(image[2], cmap=\"Blues\")\n    subax[3].set_title(\"Blue Channel\")\n    subax[4].imshow(image[3], cmap=\"Oranges\")\n    subax[4].set_title(\"Yellow Channel\")\n    return subax","metadata":{"execution":{"iopub.status.busy":"2024-10-14T11:32:12.130561Z","iopub.execute_input":"2024-10-14T11:32:12.130888Z","iopub.status.idle":"2024-10-14T11:32:12.141003Z","shell.execute_reply.started":"2024-10-14T11:32:12.130861Z","shell.execute_reply":"2024-10-14T11:32:12.140176Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"N_SAMPLES = 3\nsamples = np.array(train_df.sample(N_SAMPLES)['Id'])\n\nfig, ax = plt.subplots(N_SAMPLES,5,figsize=(20,5*N_SAMPLES))\nif ax.shape == (N_SAMPLES,):\n    ax = ax.reshape(1,-1)\nfor n in range(N_SAMPLES):\n    make_image_row(open_rgby(TRAIN, samples[n]), ax[n], samples[n].split('-')[0])","metadata":{"execution":{"iopub.status.busy":"2024-10-14T11:32:12.142072Z","iopub.execute_input":"2024-10-14T11:32:12.142384Z","iopub.status.idle":"2024-10-14T11:32:14.965094Z","shell.execute_reply.started":"2024-10-14T11:32:12.142358Z","shell.execute_reply":"2024-10-14T11:32:14.964171Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<a id='step3'></a>\n# 3. 构建PyTorch数据集\n\n在定义图片样本的读取函数之后，便可以开始构建用于PyTorch模型训练的数据集接口。由于这里的图片数据并不是标准的RGB格式，因此我们需要自定义`Dataset`，详细的构建方法可以参见此[博文](https://blog.csdn.net/u012436149/article/details/69061711)。标准的构建格式如下：\n\n```python\nclass CustomDataset(Dataset): # from torch.utils.data import Dataset\n    def __init__(self):\n        # TODO\n        # 1. Initialize file path or list of file names.\n        pass\n    def __getitem__(self, index):\n        # TODO\n        # 1. Read one data from file (e.g. using numpy.fromfile, PIL.Image.open).\n        # 2. Preprocess the data (e.g. torchvision.Transform).\n        # 3. Return a data pair (e.g. image and label).\n        pass\n    def __len__(self):\n        # You should change 0 to the total size of your dataset.\n        return 0\n```\n其中在`__init__`中完成类的继承和变量的定义，在`__getitem__`中完成对于某个index的样本的处理，并且返回该样本及其标签，在`__len__`中返回样本总量。在此，我们利用之前定义好的`train_df`、`valid_df`和`sub_df`完成数据集的构建。","metadata":{}},{"cell_type":"code","source":"class AtlasDataset(Dataset):\n    def __init__(self, df, path, size=None, label=True):        \n        self.df = df.copy()\n        self.path = path\n        self.size = size\n        self.label = label\n        if self.label:\n            self.df['Target'] = [[int(i) for i in s.split()] for s in self.df['Target']] # 对标签进行预处理\n        \n    def __getitem__(self, index):        \n        img = open_rgby(self.path, self.df['Id'].iloc[index], self.size)\n        if self.label:\n            target = np.eye(N_CLASSES,dtype=np.float)[self.df['Target'].iloc[index]].sum(axis=0) # 对标签进行独热编码\n        else:\n            target = np.zeros(N_CLASSES,dtype=np.int)\n        return img, target\n    \n    def __len__(self):\n        return len(self.df)","metadata":{"execution":{"iopub.status.busy":"2024-10-14T11:32:14.966387Z","iopub.execute_input":"2024-10-14T11:32:14.966703Z","iopub.status.idle":"2024-10-14T11:32:14.976554Z","shell.execute_reply.started":"2024-10-14T11:32:14.966674Z","shell.execute_reply":"2024-10-14T11:32:14.975439Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"在定义好`Dataset`之后，便可以使用`torch.utils.data.DataLoader`完成对数据集的封装。","metadata":{}},{"cell_type":"code","source":"size= 256\nbs = 64 # batch_size","metadata":{"execution":{"iopub.status.busy":"2024-10-14T11:32:14.977910Z","iopub.execute_input":"2024-10-14T11:32:14.978299Z","iopub.status.idle":"2024-10-14T11:32:14.988489Z","shell.execute_reply.started":"2024-10-14T11:32:14.978261Z","shell.execute_reply":"2024-10-14T11:32:14.987482Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_loader = DataLoader(AtlasDataset(train_df, TRAIN, size), batch_size=bs, shuffle=True, num_workers=2)\nvalid_loader = DataLoader(AtlasDataset(valid_df, TRAIN, size), batch_size=bs, shuffle=True, num_workers=2)\ntest_loader  = DataLoader(AtlasDataset(sub_df, TEST, size, False), batch_size=bs, shuffle=False, num_workers=2)","metadata":{"execution":{"iopub.status.busy":"2024-10-14T11:32:14.990013Z","iopub.execute_input":"2024-10-14T11:32:14.990431Z","iopub.status.idle":"2024-10-14T11:32:15.060007Z","shell.execute_reply.started":"2024-10-14T11:32:14.990391Z","shell.execute_reply":"2024-10-14T11:32:15.059177Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<a id='step4'></a>\n# 4. 搭建四通道CNN模型\n\n传统的预训练CNN模型都是基于标准的RGB三通道格式的图片数据训练得到的。然而，此竞赛的数据为四通道图片，因此不能简单调用预训练模型进行训练。为了解决这一问题，我们可以将输入层卷积层的卷积核扩增一个维度，从而达到处理四通道图片数据的目的。\n\n其中，RGB三个通道对应的卷积核参数仍采用imagenet预训练的权值，Y通道则采用随机权值。其他层仍采用imagenet权值进行初始化。\n\n## CNN模型搭建\n\n类似于数据集的搭建，PyTorch也是采用类的方法进行神经网络的搭建，其标准格式如下：\n```python\nclass CustomNet(nn.Module):\n    def __init__(self):\n        pass\n    \n    def forward(self, x):\n        return output\n\n```\n其中在`__init__`中完成类的继承和各个网络层的定义，在`forward`中完成对输入`x`前向传播的过程。PyTorch中的网络层相关接口可以查阅PyTorch官方文档：[https://pytorch.org/docs/stable/nn.html](https://pytorch.org/docs/stable/nn.html)（多看文档，熟能生巧）。","metadata":{}},{"cell_type":"code","source":"class AvgPool(nn.Module):\n    def forward(self, x):\n        return torch.squeeze(F.avg_pool2d(x, x.shape[2:]))\n\nclass ResNet50(nn.Module):\n    def __init__(self, num_classes, pretrained=True):\n        super().__init__()\n        encoder = resnet50(pretrained=False)\n        if pretrained:            \n            encoder.load_state_dict(torch.load('../input/pytorch-pretrained-models/resnet50-19c8e357.pth'))\n        \n        self.conv1 = nn.Conv2d(4, 64, kernel_size=7, stride=2, padding=3, bias=False) # 四通道卷积核\n        \n        if (pretrained):\n            w = encoder.conv1.weight\n            self.conv1.weight = nn.Parameter(torch.cat((w,0.5*(w[:,:1,:,:]+w[:,2:,:,:])),dim=1)) # 增加Y通道权值\n        \n        self.bn1 = encoder.bn1\n        self.relu = nn.ReLU(inplace=True) \n        self.maxpool = nn.MaxPool2d(kernel_size=3, stride=2, padding=1)\n        self.layer0 = nn.Sequential(self.conv1,self.relu,self.bn1,self.maxpool)\n        self.layer1 = encoder.layer1\n        self.layer2 = encoder.layer2\n        self.layer3 = encoder.layer3\n        self.layer4 = encoder.layer4\n        self.avgpool = AvgPool()\n        self.fc = nn.Sequential(nn.Dropout(p=0.5), nn.Linear(encoder.fc.in_features, num_classes)) # 引入全连接层，完成分类\n        \n    def forward(self, x):\n        x = self.layer0(x)\n        x = self.layer1(x)\n        x = self.layer2(x)\n        x = self.layer3(x)\n        x = self.layer4(x)\n        x = self.avgpool(x)\n        x = self.fc(x)        \n        return x","metadata":{"execution":{"iopub.status.busy":"2024-10-14T11:32:15.061203Z","iopub.execute_input":"2024-10-14T11:32:15.061487Z","iopub.status.idle":"2024-10-14T11:32:15.075905Z","shell.execute_reply.started":"2024-10-14T11:32:15.061460Z","shell.execute_reply":"2024-10-14T11:32:15.075040Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model = ResNet50(num_classes=N_CLASSES)","metadata":{"execution":{"iopub.status.busy":"2024-10-14T11:32:15.077217Z","iopub.execute_input":"2024-10-14T11:32:15.077571Z","iopub.status.idle":"2024-10-14T11:32:16.879333Z","shell.execute_reply.started":"2024-10-14T11:32:15.077539Z","shell.execute_reply":"2024-10-14T11:32:16.878456Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<a id='step5'></a>\n# 5. Focal Loss\n\n注意到此竞赛的数据集是一个非均衡数据集，那么如果使用一般的`BCE Loss`，由于数据类别分布原因，模型会在预测时自动对数据量大的类别产生偏移。为了解决数据分布不均导致的问题，Facebook人工智能科学家Kaiming He提出一种新的损失函数：Focal Loss，这个损失函数是在标准交叉熵损失基础上修改得到的。这个函数可以通过减少易分类样本的权重：\n\n${\\rm FL}(p_t)=-(1-p_t)^\\gamma\\log(p_t)$\n\n使得模型在训练时更专注于难分类的样本。对于Focal Loss的原理可以参考原论文（[链接](https://arxiv.org/abs/1708.02002)）。这里我们直接给出Focal Loss的PyTorch形式，其中 $\\gamma$ 取 2。","metadata":{}},{"cell_type":"code","source":"class FocalLoss(nn.Module):\n    def __init__(self, gamma=2):\n        super().__init__()\n        self.gamma = gamma\n\n    def forward(self, logit, target):\n        target = target.float()\n        max_val = (-logit).clamp(min=0)\n        loss = logit - logit * target + max_val + \\\n               ((-max_val).exp() + (-logit - max_val).exp()).log()\n\n        invprobs = F.logsigmoid(-logit * (target * 2.0 - 1.0))\n        loss = (invprobs * self.gamma).exp() * loss\n        if len(loss.size())==2:\n            loss = loss.sum(dim=1)\n        return loss.mean()","metadata":{"execution":{"iopub.status.busy":"2024-10-14T11:32:16.880541Z","iopub.execute_input":"2024-10-14T11:32:16.880868Z","iopub.status.idle":"2024-10-14T11:32:16.889154Z","shell.execute_reply.started":"2024-10-14T11:32:16.880838Z","shell.execute_reply":"2024-10-14T11:32:16.887961Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<a id='question3'></a>\n### __问题 3:__\n\n简要概述Focal Loss与BCE Loss的差别。\n\n__回答:__ ","metadata":{}},{"cell_type":"markdown","source":"Focal Loss相比于BCE Loss增加了每一项之前的权重因子，通过这种方式降低出现概率高的样本类别对于损失的贡献，从而避免类别失衡的问题。","metadata":{}},{"cell_type":"markdown","source":"<a id='step6'></a>\n# 6. 模型训练\n\n## 6.1 准备工作","metadata":{}},{"cell_type":"code","source":"os.makedirs('models') # 用于储存模型\ndevice = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\nprint(device)","metadata":{"execution":{"iopub.status.busy":"2024-10-14T11:32:16.890734Z","iopub.execute_input":"2024-10-14T11:32:16.891251Z","iopub.status.idle":"2024-10-14T11:32:16.962850Z","shell.execute_reply.started":"2024-10-14T11:32:16.891164Z","shell.execute_reply":"2024-10-14T11:32:16.961821Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"criterion = FocalLoss()\nmodel = model.to(device)\n\nlr = 0.001\noptimizer = optim.Adam(model.parameters(), lr=lr)","metadata":{"execution":{"iopub.status.busy":"2024-10-14T11:32:16.964175Z","iopub.execute_input":"2024-10-14T11:32:16.964506Z","iopub.status.idle":"2024-10-14T11:32:21.462046Z","shell.execute_reply.started":"2024-10-14T11:32:16.964476Z","shell.execute_reply":"2024-10-14T11:32:21.461100Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"PyTorch模型训练的标准格式为：\n```python\nmodel.train() # 开启训练模式，此时模型的参数可以修改\nfor batch_x, batch_y in DataLoader: # 从数据集中抽取batch\n    optimizer.zero_grad() # 清空梯度\n    output = model(batch_x) # 前向传播\n    loss = criterion(output, batch_y) # 计算损失函数\n    loss.backward() # 反向传播计算梯度\n    optimizer.step() # 更新参数\n```\n灵活应用`optimizer`的梯度操作，可以达到梯度累积`gradient accumulation`的效果。在数据较大，GPU显存不够时，可以变相增大batch_size。","metadata":{}},{"cell_type":"code","source":"def train_model(epoch, history=None):\n    model.train() \n    t = tqdm(train_loader)\n    \n    for batch_idx, (img_batch, label_batch) in enumerate(t):\n        img_batch = img_batch.to(device)\n        label_batch = label_batch.to(device)\n        \n        optimizer.zero_grad()\n        output = model(img_batch)\n        loss = criterion(output, label_batch)\n        t.set_description(f'train_loss (l={loss:.4f})')\n        \n        if history is not None:\n            history.loc[epoch + batch_idx / len(train_loader), 'train_loss'] = loss.data.cpu().numpy()\n        \n        loss.backward()    \n        optimizer.step()\n    \n    torch.save(model.state_dict(), 'models/epoch{}.pth'.format(epoch))","metadata":{"execution":{"iopub.status.busy":"2024-10-14T11:32:21.463423Z","iopub.execute_input":"2024-10-14T11:32:21.463851Z","iopub.status.idle":"2024-10-14T11:32:21.471897Z","shell.execute_reply.started":"2024-10-14T11:32:21.463807Z","shell.execute_reply":"2024-10-14T11:32:21.470663Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def binarize_prediction(probabilities, threshold: float, argsorted=None,\n                        min_labels=1, max_labels=10):\n    \"\"\" Return matrix of 0/1 predictions, same shape as probabilities.\n    \"\"\"\n    assert probabilities.shape[1] == N_CLASSES\n    if argsorted is None:\n        argsorted = probabilities.argsort(axis=1)\n    max_mask = _make_mask(argsorted, max_labels)\n    min_mask = _make_mask(argsorted, min_labels)\n    prob_mask = probabilities > threshold\n    return (max_mask & prob_mask) | min_mask\n\ndef _make_mask(argsorted, top_n: int):\n    mask = np.zeros_like(argsorted, dtype=np.uint8)\n    col_indices = argsorted[:, -top_n:].reshape(-1)\n    row_indices = [i // top_n for i in range(len(col_indices))]\n    mask[row_indices, col_indices] = 1\n    return mask","metadata":{"execution":{"iopub.status.busy":"2024-10-14T11:32:21.473519Z","iopub.execute_input":"2024-10-14T11:32:21.474010Z","iopub.status.idle":"2024-10-14T11:32:21.483964Z","shell.execute_reply.started":"2024-10-14T11:32:21.473969Z","shell.execute_reply":"2024-10-14T11:32:21.482922Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"我们需要设定阈值，使得模型可以根据输出值与阈值的相对大小判断图像中是否存在某类蛋白。显然该阈值无法人为确定，这里我们可以观察不同阈值在验证集上的表现（Macro F1），从而确定最佳的阈值。","metadata":{}},{"cell_type":"code","source":"def evaluate(epoch, history=None): # 验证函数\n    model.eval() # 开启验证模式，此时模型的参数不可修改\n    valid_loss = 0.\n    all_predictions, all_targets = [], []\n    \n    with torch.no_grad():\n        for batch_idx, (img_batch, label_batch) in enumerate(valid_loader):\n            all_targets.append(label_batch.numpy().copy())\n            img_batch = img_batch.to(device)\n            label_batch = label_batch.to(device)\n\n            output = model(img_batch)\n            loss = criterion(output, label_batch)\n            valid_loss += loss.data\n            predictions = torch.sigmoid(output)\n            all_predictions.append(predictions.cpu().numpy())\n    all_predictions = np.concatenate(all_predictions)\n    all_targets = np.concatenate(all_targets)\n    \n    valid_loss /= (batch_idx+1)\n    \n    if history is not None:\n        history.loc[epoch, 'valid_loss'] = valid_loss.cpu().numpy()\n    \n    print('Epoch: {}\\tLR: {:.6f}\\tValid Loss: {:.4f}'.format(\n        epoch, optimizer.state_dict()['param_groups'][0]['lr'], valid_loss))\n    \n    def get_score(y_pred):\n        with warnings.catch_warnings():\n            warnings.simplefilter('ignore', category=UndefinedMetricWarning)\n            return f1_score(all_targets, y_pred, average='macro')\n    \n    metrics = {}\n    argsorted = all_predictions.argsort(axis=1)\n    # 测试不同阈值\n    for threshold in [0.05,0.10,0.15,0.2,0.25,0.3,0.35,0.4,0.45]: \n        metrics[threshold] = get_score(\n            binarize_prediction(all_predictions, threshold, argsorted))\n    best_thr = max(metrics, key=metrics.get)\n    print(' | '.join(f'thr_{k:.2f} {v:.3f}' for k, v in sorted(\n        metrics.items(), key=lambda kv: -kv[1])[:5]))\n    \n    return valid_loss, best_thr, metrics[best_thr]","metadata":{"execution":{"iopub.status.busy":"2024-10-14T11:32:21.485664Z","iopub.execute_input":"2024-10-14T11:32:21.486074Z","iopub.status.idle":"2024-10-14T11:32:21.501406Z","shell.execute_reply.started":"2024-10-14T11:32:21.486034Z","shell.execute_reply":"2024-10-14T11:32:21.500531Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 6.2 开始训练","metadata":{}},{"cell_type":"code","source":"history_train = pd.DataFrame()\nhistory_valid = pd.DataFrame()\n\nn_epochs = 100\ninit_epoch = 0\nmax_lr_changes = 1\nvalid_losses = []\nthreshold = {}\nmacro_f1 = {}\nlr_reset_epoch = init_epoch\npatience = 2\nlr_changes = 0\nbest_valid_loss = 1000.\n\nfor epoch in range(init_epoch, n_epochs):\n    torch.cuda.empty_cache()\n    gc.collect()\n    train_model(epoch, history_train)\n    valid_loss, best_thr, best_f1 = evaluate(epoch, history_valid)\n    valid_losses.append(valid_loss)\n    threshold[epoch] = best_thr\n    macro_f1[epoch] = best_f1\n\n    if valid_loss < best_valid_loss:\n        best_valid_loss = valid_loss\n    elif (patience and epoch - lr_reset_epoch > patience and\n          min(valid_losses[-patience:]) > best_valid_loss):\n        # \"patience\" epochs without improvement\n        lr_changes +=1\n        if lr_changes > max_lr_changes: # 早期停止\n            break\n        lr /= 5 # 学习率衰减\n        print(f'lr updated to {lr}')\n        lr_reset_epoch = epoch\n        optimizer.param_groups[0]['lr'] = lr","metadata":{"execution":{"iopub.status.busy":"2024-10-14T11:32:21.505089Z","iopub.execute_input":"2024-10-14T11:32:21.505385Z","iopub.status.idle":"2024-10-14T13:05:40.175589Z","shell.execute_reply.started":"2024-10-14T11:32:21.505357Z","shell.execute_reply":"2024-10-14T13:05:40.174363Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 6.3 训练过程可视化","metadata":{}},{"cell_type":"code","source":"fig, ax = plt.subplots(1,3,figsize=(20,5))\n\nax[0].plot(history_train['train_loss'].iloc[100:])\nax[0].set_xlabel('Epoch')\nax[0].set_ylabel('Train Loss')\n\nax[1].plot(history_valid['valid_loss'])\nax[1].set_xlabel('Epoch')\nax[1].set_ylabel('Valid Loss')\n\nax[2].scatter(macro_f1.keys(), macro_f1.values())\nax[2].set_xlabel('Epoch')\nax[2].set_ylabel('Macro F1 Score')","metadata":{"execution":{"iopub.status.busy":"2024-10-14T13:05:40.177764Z","iopub.execute_input":"2024-10-14T13:05:40.178055Z","iopub.status.idle":"2024-10-14T13:05:40.653124Z","shell.execute_reply.started":"2024-10-14T13:05:40.178023Z","shell.execute_reply":"2024-10-14T13:05:40.652186Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<a id='step7'></a>\n# 7. 后处理与提交","metadata":{}},{"cell_type":"code","source":"# 加载最佳模型\nbest_epoch = max(macro_f1, key=macro_f1.get)\nprint('The best epoch is epoch {}'.format(best_epoch))\n\nmodel.load_state_dict(torch.load('models/epoch{}.pth'.format(best_epoch)))\nmodel.eval();","metadata":{"execution":{"iopub.status.busy":"2024-10-14T13:05:40.654402Z","iopub.execute_input":"2024-10-14T13:05:40.654702Z","iopub.status.idle":"2024-10-14T13:05:40.769657Z","shell.execute_reply.started":"2024-10-14T13:05:40.654672Z","shell.execute_reply":"2024-10-14T13:05:40.768579Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Inference\noutputlist = []\nfor img_batch, _ in tqdm(test_loader):\n    with torch.no_grad():\n        output = torch.sigmoid(model(img_batch.to(device)))\n    output = output.data.cpu().numpy()\n    for i in output: \n        outputlist.append(i)","metadata":{"execution":{"iopub.status.busy":"2024-10-14T13:05:40.771083Z","iopub.execute_input":"2024-10-14T13:05:40.771488Z","iopub.status.idle":"2024-10-14T13:11:02.726498Z","shell.execute_reply.started":"2024-10-14T13:05:40.771444Z","shell.execute_reply":"2024-10-14T13:11:02.725354Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"thr = threshold[best_epoch]\nprint('The best threshold is {:.3f}'.format(thr))\n\nprediction = [' '.join([str(i) for i in np.argwhere((j > thr).astype(int)==1).reshape(-1)]) for j in outputlist]\nsub_df['Predicted'] = prediction\nsub_df.to_csv('submission.csv', index=False)","metadata":{"execution":{"iopub.status.busy":"2024-10-14T13:11:02.728204Z","iopub.execute_input":"2024-10-14T13:11:02.728550Z","iopub.status.idle":"2024-10-14T13:11:03.385435Z","shell.execute_reply.started":"2024-10-14T13:11:02.728515Z","shell.execute_reply":"2024-10-14T13:11:03.384548Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}