{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport os, glob, pickle, time, gc, copy, sys\nimport warnings\nimport requests\n\nimport pydicom\nimport cv2\nimport os, os.path as osp\n\nwarnings.filterwarnings('ignore')\npd.set_option('display.max_columns', 100)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"# load train data\ndf_train = pd.read_csv(\"../input/rsna-str-pulmonary-embolism-detection/train.csv\")\nprint(df_train.shape)\ndf_train.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# extract exam (study) level data\ncol_index = 'SOPInstanceUID'\ncol_groupby = 'StudyInstanceUID'\ndf_train_study = df_train[df_train[col_groupby].duplicated()==False].reset_index(drop=True)\nprint(df_train_study.shape)\ndf_train_study.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# calculate q_i\ndf_tmp = df_train.groupby(col_groupby)['pe_present_on_image'].agg(len).reset_index()\ndf_tmp.columns = [col_groupby, 'num_images']\ndf_train_study = pd.merge(df_train_study, df_tmp, on=col_groupby, how='left')\n\ndf_tmp = df_train.groupby(col_groupby)['pe_present_on_image'].agg('sum').reset_index()\ndf_tmp.columns = [col_groupby, 'm_i']\ndf_train_study = pd.merge(df_train_study, df_tmp, on=col_groupby, how='left')\n\ndf_train_study['q_i'] = df_train_study['m_i'] /df_train_study['num_images'] \ndf_train = pd.merge(df_train, df_train_study[[col_groupby, 'num_images', 'q_i']], on=col_groupby, how='left')\ndf_train.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Mean prediction\nThis competition's metric is logloss. \nThe representative value minimizes logloss is average. "},{"metadata":{"trusted":true},"cell_type":"code","source":"# calculate average\ncol_index = 'SOPInstanceUID'\ncol_targets = [\n    'negative_exam_for_pe',\n    'indeterminate',\n    'chronic_pe',\n    'acute_and_chronic_pe',\n    'central_pe',\n    'leftsided_pe',\n    'rightsided_pe',\n    'rv_lv_ratio_gte_1',\n    'rv_lv_ratio_lt_1',\n    'pe_present_on_image',\n]\nmean_targets = np.zeros(len(col_targets), np.float32)\nfor i, col in enumerate(col_targets[:-1]):\n    mean_targets[i] = df_train_study[col].mean()\nmean_targets[-1] = df_train[col_targets[-1]].mean()\npreds_mean_study = np.ones([len(df_train_study), len(col_targets[:-1])], np.float32) * mean_targets[:-1][np.newaxis]\npreds_mean_image = np.ones(len(df_train), np.float32) * mean_targets[-1]\n\nfor i, col in enumerate(col_targets):\n    print(\"{} average: {:.6f}\".format((col +\" \"*50)[:30], mean_targets[i]))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# calculate metrics\nfrom sklearn import metrics\ndef calc_metrics(y_true_exam, y_pred_exam, y_true_imag, y_pred_imag, q_image):\n    weights = np.array([\n        0.0736196319, \n        0.09202453988, \n        0.1042944785, \n        0.1042944785, \n        0.1877300613, \n        0.06257668712, \n        0.06257668712,\n        0.2346625767,\n        0.0782208589,\n        0.07361963,\n    ])\n    score_list = []\n    scores = {}\n    for i in range(9):\n        bce = metrics.log_loss(y_true_exam[:,i], y_pred_exam[:,i])\n        scores[col_targets[i]] = bce\n        score_list.append(bce)\n        print(\"{}: {:.6f}\".format((col_targets[i]+\" \"*50)[:30], bce))\n    score_s = np.sum(weights[:-1]*np.array(score_list))\n    \n    scores[\"exam_level_weighted_log_loss\"] = score_s\n    \n    print(\"{}: {:.6f}\".format((\"exam_level_weighted_log_loss\" +\" \"*50)[:30], score_s))\n    score_i =  np.sum(- q_image * (y_true_imag*np.log(y_pred_imag) + (1-y_true_imag)*np.log(1-y_pred_imag))) / np.sum(q_image)\n    scores[\"image_level_weighted_log_loss\"] = score_i\n    print(\"{}: {:.6f}\".format((\"image_level_weighted_log_loss\" +\" \"*50)[:30], score_i))\n#     scores.append(score)\n    score_all = (score_s * len(y_true_exam) + score_i * np.sum(q_image) * weights[-1]) / (len(y_true_exam)+np.sum(q_image)* weights[-1])\n    scores[\"total_loss\"] = score_all\n    print(\"{}: {:.6f}\".format((\"total_loss\" +\" \"*50)[:30], score_all))\n#     print(len(y_true_exam), np.sum(q_image)* weights[-1])\n    return scores\n\n\n_ = calc_metrics(\n    y_true_exam=df_train_study[col_targets[:-1]].values, \n    y_pred_exam=preds_mean_study, \n    y_true_imag=df_train[col_targets[-1]].values, \n    y_pred_imag=preds_mean_image, \n    q_image=df_train['q_i'].values)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This prediction get a score of 0.60 on the train data and 0.56 on the pulic leaderboard."},{"metadata":{},"cell_type":"markdown","source":"# Weighted mean prediction\nActually, the metric for image-level prediction is **weighted** logloss. \nThe representative value minimizes weighted logloss is **weighted** average."},{"metadata":{"trusted":true},"cell_type":"code","source":"q_weighted_mean = np.sum(df_train['pe_present_on_image'] * df_train['q_i']) / np.sum(df_train['q_i'])\nprint(\"q_weighted_mean: {:.6f}\".format(q_weighted_mean))\npreds_mean_image_q_weighted = np.ones(len(df_train), np.float32) * q_weighted_mean","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"_ = calc_metrics(\n    y_true_exam=df_train_study[col_targets[:-1]].values, \n    y_pred_exam=preds_mean_study, \n    y_true_imag=df_train[col_targets[-1]].values, \n    y_pred_imag=preds_mean_image_q_weighted, \n    q_image=df_train['q_i'].values)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This prediction get a score of 0.46 on the train data and 0.44 on the pulic leaderboard."},{"metadata":{},"cell_type":"markdown","source":"# Weighted mean prediction per slice location\nWe can get farther. \npulmonary embolism is located on pulmonary artery ,so that pe_present_on_image must be correlated with slice location."},{"metadata":{"trusted":true},"cell_type":"code","source":"# get dicom paths\ndf_train['path'] = (\"../input/rsna-str-pulmonary-embolism-detection/train/\" \n                   + df_train['StudyInstanceUID'].values + \"/\"\n                   + df_train['SeriesInstanceUID'].values + \"/\"\n                   + df_train['SOPInstanceUID'].values + \".dcm\"\n                  )\nprint(df_train['path'][0])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# get series index of image\nimport multiprocessing\nfrom concurrent.futures import ProcessPoolExecutor\n\ndef task(i):\n    if (i+1)%100000==0:\n        print(\"{}/{} {:.1f}\".format(i+1, len(df_train), time.time()-starttime))\n    path = df_train['path'][i]\n    tmp_dcm = pydicom.dcmread(path)\n    return tmp_dcm.ImagePositionPatient[-1]\n\n\nstarttime = time.time()\nexecutor = ProcessPoolExecutor(max_workers=multiprocessing.cpu_count())\n# futures = [executor.submit(task, i) for i in range(10000)]\nfutures = [executor.submit(task, i) for i in range(len(df_train.iloc[:]))]\nresult_list = []\nfor i in range(len(futures)):\n    result_list.append(futures[i].result())\ndf_train['z_pos'] = 0\ndf_train['z_pos'][:len(result_list)] = result_list\n\ndel futures, result_list\ngc.collect()\n\ndf_train.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# calculate slice location\ndf_tmp = []\nfor i in range(len(df_train_study)):\n    if (i+1)%1000==0: print(\"{}/{}\".format(i+1, len(df_train_study)))\n    study = df_train_study[col_groupby][i]\n    df_study = df_train[df_train[col_groupby]==study].sort_values('z_pos').reset_index(drop=True)\n    df_study['series_index'] = np.arange(len(df_study))\n    df_tmp.append(df_study[[col_index, 'series_index']])\ndf_tmp = pd.concat(df_tmp)\n\ndf_train = pd.merge(df_train, df_tmp, on=col_index, how='left')\n# df_test = pd.merge(df_test, df_test_study[[col_groupby, 'num_images']], on=col_groupby, how='left')\ndf_train['slice_location'] = df_train['series_index'] / (df_train['num_images'] - 1)\ndf_train.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# visualize the relation between pe and slice location\nplt.figure(figsize=(10, 5))\nplt.hist(df_train['slice_location'][df_train['pe_present_on_image']==True], bins=np.arange(101)/100, label='PE', density=True, alpha=0.3)\nplt.hist(df_train['slice_location'][df_train['pe_present_on_image']==False], bins=np.arange(101)/100, label='no PE', density=True, alpha=0.3)\nplt.legend()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Obviously, there is a strong correlation between PE and slice location. \nTherefore, the weighted average per slice location must be a good representative value."},{"metadata":{"trusted":true},"cell_type":"code","source":"bins = 8\ndf_train['bins'] = bins-1\nfor i in range(bins):\n    df_train['bins'][(df_train['slice_location']>=(i/bins)) & (df_train['slice_location']<((i+1)/bins))] = i\ndf_train.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_train['slice_location'].min(), df_train['slice_location'].max()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for i in range(bins):\n    print(i, np.sum(df_train['bins']==i))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"q_weighted_means = np.zeros(bins, np.float32)\nfor i in range(bins):\n    tmp_index = df_train['bins']==i\n    q_weighted_means[i] = np.sum(df_train['pe_present_on_image'][tmp_index].values * df_train['q_i'][tmp_index].values) / np.sum(df_train['q_i'][tmp_index].values)\ndf_train['q_weighted_means'] = df_train['bins'].apply(lambda x: q_weighted_means[x])\nprint(q_weighted_means)\ndf_train.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"q_weighted_means = np.array([0.00326324, 0.05970682, 0.32645303, 0.67452216, 0.71344817, 0.4734337, 0.0740926, 0.00369781])\nprint(q_weighted_means)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"_ = calc_metrics(\n    y_true_exam=df_train_study[col_targets[:-1]].values, \n    y_pred_exam=preds_mean_study, \n    y_true_imag=df_train[col_targets[-1]].values, \n    y_pred_imag=df_train['q_weighted_means'].values, \n    q_image=df_train['q_i'].values)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This prediction get a score of 0.35 on the train data and 0.33 on the pulic leaderboard."},{"metadata":{},"cell_type":"markdown","source":"# Make submission"},{"metadata":{"trusted":true},"cell_type":"code","source":"# load test data\ndf_test = pd.read_csv(\"../input/rsna-str-pulmonary-embolism-detection/test.csv\")\nprint(df_test.shape)\ndf_test.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# get dicom paths\ndf_test['path'] = (\"../input/rsna-str-pulmonary-embolism-detection/test/\" \n                   + df_test['StudyInstanceUID'].values + \"/\"\n                   + df_test['SeriesInstanceUID'].values + \"/\"\n                   + df_test['SOPInstanceUID'].values + \".dcm\"\n                  )\nprint(df_test['path'][0])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# extract exam (study) level data\ndf_test_study = df_test[df_test[col_groupby].duplicated()==False].reset_index(drop=True)\ndf_tmp = df_test.groupby(col_groupby)[col_index].agg(len).reset_index()\ndf_tmp.columns = [col_groupby, 'num_images']\ndf_test_study = pd.merge(df_test_study, df_tmp, on=col_groupby, how='left')\ndf_test = pd.merge(df_test, df_test_study[[col_groupby, 'num_images']], on=col_groupby, how='left')\nprint(df_test.shape)\ndf_test.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# get series index of image\ndef task(i):\n    if (i+1)%10000==0:\n        print(\"{}/{} {:.1f}\".format(i+1, len(df_test), time.time()-starttime))\n    path = df_test['path'][i]\n    tmp_dcm = pydicom.dcmread(path)\n    return tmp_dcm.ImagePositionPatient[-1]\n\nstarttime = time.time()\nexecutor = ProcessPoolExecutor(max_workers=multiprocessing.cpu_count())\n# futures = [executor.submit(task, i) for i in range(10000)]\nfutures = [executor.submit(task, i) for i in range(len(df_test))]\nresult_list = []\nfor i in range(len(futures)):\n    result_list.append(futures[i].result())\ndf_test['z_pos'] = result_list\ndf_test.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# calculate slice location\ndf_tmp = []\nfor i in range(len(df_test_study)):\n    if (i+1)%100==0: print(\"{}/{}\".format(i+1, len(df_test_study)))\n    study = df_test_study[col_groupby][i]\n    df_study = df_test[df_test[col_groupby]==study].sort_values('z_pos').reset_index(drop=True)\n    df_study['series_index'] = np.arange(len(df_study))\n    df_tmp.append(df_study[[col_index, 'series_index']])\ndf_tmp = pd.concat(df_tmp)\n\ndf_test = pd.merge(df_test, df_tmp, on=col_index, how='left')\n# df_test = pd.merge(df_test, df_test_study[[col_groupby, 'num_images']], on=col_groupby, how='left')\ndf_test['slice_location'] = df_test['series_index'] / (df_test['num_images'] - 1)\ndf_test.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# get weighted mean prediction per slice location\nbins = 8\ndf_test['bins'] = bins-1\nfor i in range(bins):\n    df_test['bins'][(df_test['slice_location']>=(i/bins)) & (df_test['slice_location']<((i+1)/bins))] = i\ndf_test['q_weighted_means'] = df_test['bins'].apply(lambda x: q_weighted_means[x])\ndf_test.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_sub_tmp = copy.deepcopy(df_test[[col_index, 'q_weighted_means']])\ndf_sub_tmp.columns = ['id', 'label']\nfor i, col in enumerate(col_targets[:-1]):\n    df_tmp = df_test_study[[col_groupby]]\n    df_tmp.columns = ['id']\n    df_tmp['label'] = mean_targets[i]\n    df_tmp['id'] = df_tmp['id'] + '_{}'.format(col)\n    df_sub_tmp = pd.concat([df_sub_tmp, df_tmp])\ndf_sub_tmp = df_sub_tmp.reset_index(drop=True)\nprint(df_sub_tmp.shape)\ndf_sub_tmp.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_sub = pd.read_csv(\"../input/rsna-str-pulmonary-embolism-detection/sample_submission.csv\")\nprint(df_sub.shape)\ndf_sub = pd.merge(df_sub[['id']], df_sub_tmp, on='id', how='left')\ndf_sub = df_sub.fillna(0.5)\ndf_sub.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_sub.to_csv(\"submission.csv\", index=None)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for i in range(bins):\n    print(i, np.sum(df_test['bins']==i))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}