{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Install dependencies"},{"metadata":{"trusted":true},"cell_type":"code","source":"'''\n!conda install -y pillow\n!conda install -y scikit-learn\n!conda install -y scikit-image\n!conda install -y tqdm\n!conda install -c conda-forge -y gdcm\n!pip install dicom-numpy\n!pip install pydicom\n'''","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Imports"},{"metadata":{"trusted":true},"cell_type":"code","source":"import os\nfrom functools import partial\nfrom functools import reduce\n\nimport numpy as np\nimport pandas as pd\nfrom tqdm import tqdm\nimport tensorflow as tf\n\nphysical_devices = tf.config.list_physical_devices('GPU')\nfor gpu_instance in physical_devices:\n    tf.config.experimental.set_memory_growth(gpu_instance, True)\n    \nfrom skimage.transform import resize\nimport pydicom","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Constants"},{"metadata":{"trusted":true},"cell_type":"code","source":"MIN_WEEK = -12\nMAX_WEEK = 133\nMAX_FVC = 4000\n#../input/osic-pulmonary-fibrosis-progression\n#INPUT_ROOT = '../input/osic-pulmonary-fibrosis-progression'\nINPUT_ROOT = '/kaggle/input/osic-pulmonary-fibrosis-progression'\n\nTRAIN_SCANS_ROOT = os.path.join(INPUT_ROOT, 'train')\nTEST_SCANS_ROOT = os.path.join(INPUT_ROOT, 'test')\nSCAN_DEPTH = 20\nHEIGHT = 32\nWIDTH = 32\nLEARNING_RATE = 0.001\n\nTRAIN_INPUT_FILE = os.path.join(INPUT_ROOT, 'train.csv')\nTEST_INPUT_FILE = os.path.join(INPUT_ROOT, 'test.csv')\nBATCH_SIZE = 192\nSEQUENCE_LENGTH = 1\n\nSTEPS_PER_EPOCH = 1500 // BATCH_SIZE\nVALIDATION_STEPS = 300 // BATCH_SIZE\n\nEPOCHS = 50\n\nTEST_OUTPUT = 'submission.csv'","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Data utility functions"},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_patient_scan(patient_dir, n_depth=5, rows=64, columns=64):\n    patient_files = [os.path.join(e) for e in os.listdir(patient_dir)]\n    patient_files.sort(key=lambda fname: int(fname.split('.')[0]))\n    dcm_slices = [pydicom.read_file(os.path.join(patient_dir, f)) for f in patient_files]\n    # Resample slices such that the depth of the CT scan is 'n_depth'\n    slice_group = n_depth / len(patient_files)\n    slice_indexes = [int(idx / slice_group) for idx in range(n_depth)]\n    dcm_slices = [dcm_slices[i] for i in slice_indexes]\n    # Merge slices\n    shape = (rows, columns)\n    shape = (n_depth, *shape)\n    img = np.empty(shape, dtype='float32')\n    for idx, dcm in enumerate(dcm_slices):\n        # Rescale and shift in order to get accurate pixel values\n        slope = float(dcm.RescaleSlope)\n        intercept = float(dcm.RescaleIntercept)\n        resized_img = resize(dcm.pixel_array.astype('float32'), (rows, columns), anti_aliasing=True)\n        img[idx, ...] = resized_img * slope + intercept\n    return img\n\n\ndef get_dicom_data(patients_root, n_depth=5, rows=64, columns=64):\n    def gen(patients_root):\n        for patient_dir in os.listdir(patients_root):\n            patient_dir = os.path.join(patients_root, patient_dir)\n            img = get_patient_scan(patient_dir, n_depth=n_depth, rows=rows, columns=columns)\n            yield img\n\n    return tf.data.Dataset.from_generator(partial(gen, patients_root), output_types=tf.float32)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Laplace log likelihood metric and loss"},{"metadata":{"trusted":true},"cell_type":"code","source":"def laplace_log_likelihood(y_true, y_pred):\n    uncertainty_clipped = tf.maximum(y_pred[:, 1:2] * 1000.0, 70)\n    prediction = y_pred[:, :1]\n    delta = tf.minimum(tf.abs(y_true - prediction), 1000.0)\n    metric = -np.sqrt(2.0) * delta / uncertainty_clipped - tf.math.log(np.sqrt(2.0) * uncertainty_clipped)\n    return tf.reduce_mean(metric)\n\ndef laplace_log_likelihood_loss(y_true, y_pred):\n    uncertainty_clipped = tf.maximum(y_pred[:, 1:2] * 1000.0, 70)\n    prediction = y_pred[:, :1]\n    delta = tf.minimum(tf.abs(y_true - prediction), 1000.0)\n    metric = -np.sqrt(2.0) * delta / uncertainty_clipped - tf.math.log(np.sqrt(2.0) * uncertainty_clipped)\n    return -tf.reduce_mean(metric)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Build combined RNN / 3D CNN model"},{"metadata":{},"cell_type":"markdown","source":""},{"metadata":{"trusted":true},"cell_type":"code","source":"def get3dcnn_model(width, height, depth):\n    # Do we have to specify channels ?\n    inputs = tf.keras.Input(shape=(width, height, depth, 1))\n    x = tf.keras.layers.Conv3D(32, kernel_size=(5, 5, 5), activation='relu')(inputs)\n    x = tf.keras.layers.MaxPool3D()(x)\n\n    x = tf.keras.layers.Conv3D(64, kernel_size=(5, 5, 5), activation='relu')(x)\n    x = tf.keras.layers.MaxPool3D()(x)\n\n    return inputs, x\n\n\n\ndef get_combined_model(sequence_length, learning_rate, width, height, depth):\n    rnn_inputs = tf.keras.Input(shape=(sequence_length, 2))\n    x = tf.keras.layers.Masking(mask_value=-1, input_shape=(sequence_length, 1))(rnn_inputs)\n    rnn_out = 4\n    rnn_out = tf.keras.layers.GRU(rnn_out)(x)\n\n    cnn3d_inputs, cnn3d_out = get3dcnn_model(width, height, depth)\n    cnn3d_out_shape = reduce(lambda x, y: x*y, cnn3d_out.shape[1:])\n    cnn3d_out = tf.keras.layers.Reshape((cnn3d_out_shape,))(cnn3d_out)\n\n    combined_out = tf.keras.layers.concatenate([rnn_out, cnn3d_out])\n\n    prediction_output = tf.keras.layers.Dense(1)(combined_out)\n    uncertainty_output = tf.keras.layers.Dense(1, activation='sigmoid')(combined_out)\n\n    outputs = tf.keras.layers.concatenate([prediction_output, uncertainty_output])\n\n    model = tf.keras.Model(inputs=[rnn_inputs, cnn3d_inputs], outputs=outputs)\n\n    metrics = [laplace_log_likelihood]\n\n    model.compile(loss=laplace_log_likelihood_loss, optimizer=tf.optimizers.Adam(learning_rate=learning_rate), metrics=metrics)\n\n    return model","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Dataset generator"},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_combined_data(input_file, batch_size, sequence_length, scans_root, n_depth, rows, columns, max_fvc, split=0.8):\n    train_data = pd.read_csv(input_file)\n    n_features = 0\n    train_data['Weeks'] = train_data['Weeks']\n    n_features += 1\n    train_data['FVC'] /= max_fvc\n    n_features += 1\n    grouped = train_data.groupby(train_data.Patient)\n    n_data = len(train_data)\n\n    def gen():\n        for patient in tqdm(train_data['Patient'].unique()):\n            patient_df = grouped.get_group(patient)\n            FVC = patient_df['FVC'].iloc[:-1].tolist()\n            weeks = patient_df['Weeks']\n            Weeks = weeks.iloc[:-1]\n            Weeks_next = weeks.iloc[1:]\n            week_diff = (np.array(Weeks_next.tolist()) - np.array(Weeks.tolist())) / (MAX_WEEK - MIN_WEEK)\n            converted_data = {'FVC': FVC, 'week_diff': week_diff}\n            converted_df = pd.DataFrame.from_dict(converted_data)\n            indexes = sorted(list(converted_df.index))\n            patient_dir = os.path.join(scans_root, patient)\n            img = np.expand_dims(get_patient_scan(patient_dir, n_depth=n_depth, rows=rows, columns=columns), axis=-1)\n\n            for idx in indexes:\n                prev_indexes = sorted(list(range(int(idx - sequence_length), int(idx))))\n                if len(set(indexes).intersection(set(prev_indexes))):\n                    sequence = np.empty((sequence_length, n_features))\n                    for i, prev_idx in enumerate(prev_indexes):\n                        if prev_idx in converted_df['FVC'].index:\n                            sequence[i] = [converted_df['FVC'].loc[prev_idx], converted_df['week_diff'].loc[prev_idx]]\n                        else:\n                            sequence[i] = [-1, -1]\n                    yield ((sequence, img), converted_df['FVC'].loc[idx])\n\n    dataset = tf.data.Dataset.from_generator(gen, output_types=((tf.float32, tf.float32), tf.float32)).repeat(None).shuffle(n_data)\n\n    train_size = int(split * n_data)\n    train_dataset = dataset.take(train_size)\n    val_dataset = dataset.skip(train_size)\n\n    return train_dataset.batch(batch_size), val_dataset.batch(batch_size)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Training model"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_dataset, val_dataset = get_combined_data(TRAIN_INPUT_FILE, BATCH_SIZE, SEQUENCE_LENGTH, TRAIN_SCANS_ROOT, SCAN_DEPTH, HEIGHT, WIDTH, max_fvc=MAX_FVC)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model = get_combined_model(SEQUENCE_LENGTH, LEARNING_RATE, WIDTH, HEIGHT, SCAN_DEPTH)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model.fit(train_dataset, epochs=EPOCHS, validation_data=val_dataset, steps_per_epoch=STEPS_PER_EPOCH, validation_steps=VALIDATION_STEPS)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Predict"},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_input_data(FVC, img, sequence, week_diff):\n    converted_data = {'FVC': FVC, 'week_diff': [week_diff]}\n    converted_df = pd.DataFrame.from_dict(converted_data)\n    sequence[0] = [converted_df['FVC'].loc[0], converted_df['week_diff'].loc[0]]\n    data = (np.expand_dims(sequence, axis=0), np.expand_dims(img, axis=0))\n    return data","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"    test_data = pd.read_csv(TEST_INPUT_FILE)\n    n_features = 0\n    test_data['Weeks'] = test_data['Weeks']\n    n_features += 1\n    test_data['FVC'] /= MAX_FVC\n    n_features += 1\n    grouped = test_data.groupby(test_data.Patient)\n\n    prediction_data = {'Patient_Week': [], 'FVC': [], 'Confidence': []}\n    all_weeks = set(list(range(MIN_WEEK, MAX_WEEK + 1)))\n\n    for patient in tqdm(test_data['Patient'].unique()):\n        patient_df = grouped.get_group(patient)\n        FVC = patient_df['FVC'].iloc[:1].tolist()\n        measurement_week = patient_df['Weeks'].iloc[0]\n        prediction_weeks = sorted(all_weeks - set([measurement_week]))\n        week_diffs = (np.array(prediction_weeks) - np.array(measurement_week)) / (MAX_WEEK - MIN_WEEK)\n        patient_dir = os.path.join(TEST_SCANS_ROOT, patient)\n        img = np.expand_dims(get_patient_scan(patient_dir, n_depth=SCAN_DEPTH, rows=HEIGHT, columns=WIDTH), axis=-1)\n        sequence = np.empty((1, n_features))\n\n        prediction_data['Patient_Week'].append(patient + '_' + str(measurement_week))\n        prediction_data['FVC'].append(int(FVC[0] * MAX_FVC))\n        prediction_data['Confidence'].append(100)\n\n        for week, week_diff in zip(prediction_weeks, week_diffs):\n            data = get_input_data(FVC, img, sequence, week_diff)\n            prediction = model.predict([data])\n            FVC = prediction[0][0]\n            uncertainty = prediction[0][1]\n            prediction_data['Patient_Week'].append(patient + '_' + str(week))\n            prediction_data['FVC'].append(int(FVC * MAX_FVC))\n            #confidence = 1/(uncertainty+1) * 100.0\n            confidence = uncertainty\n            prediction_data['Confidence'].append(confidence)\n\n    indexes = list(range(len(prediction_data['Patient_Week'])))\n\n    def get_key(patient_week):\n        patient, week = patient_week.split('_')\n        return int(week)\n\n    sorted_data = sorted(zip(indexes, prediction_data['Patient_Week'], prediction_data['FVC'], prediction_data['Confidence']), key=lambda e: get_key(e[1]))\n    prediction_data['Patient_Week'] = [e[1] for e in sorted_data]\n    prediction_data['FVC'] = [e[2] for e in sorted_data]\n    prediction_data['Confidence'] = [e[3] for e in sorted_data]\n\n    df = pd.DataFrame.from_dict(prediction_data)\n    df.to_csv(TEST_OUTPUT, index=False)","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}