{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Imports and hyperparameters"},{"metadata":{"trusted":true},"cell_type":"code","source":"import os, pydicom, random, math\nimport pandas as pd\nimport numpy as np\nimport tensorflow as tf\nimport matplotlib.pyplot as plt\n%matplotlib inline\nfrom sklearn.preprocessing import OneHotEncoder, StandardScaler, PolynomialFeatures\nfrom sklearn.cluster import KMeans\nfrom sklearn.linear_model import LinearRegression\nfrom scipy.optimize import minimize\n\nDATA_PATH = \"../input/osic-pulmonary-fibrosis-progression/\"\nN_FOLDS = 5                                   # Number of component models/data splits\nMIN_TEST_WEEK = -12\nMAX_TEST_WEEK = 133\nTARGET_VAR = 'targetFVC'\nMAX_SMOOTHING_RANGE = 5                       # Range for observation smoothing\n\n# Image processing\nMAX_PLANES, MAX_ROW, MAX_COL = 10, 100, 100   # Maximum dimensions on Z-, X- and Y- axes\nINVERSE_WEIGHT_FCT = lambda y: y**2           # Inverse distance penalty function for weighting\nN_NEIGHBOURS = 8                              # Number of planes to consider for synthetic planes\nPIXEL_VALUE_RANGE = 32747 + 15000\n\n# Model and training hyperparameters\nno_img_params = {'shapes': [300,80]}                               # Tabular data model\nw_img_params = {'shapes': [700,400,50], 'kernels': 10}             # Tabular+imaging data model\nfit_params = {'sample_size': 5000, 'epochs': 10, 'batch_size': 10} # Training\n\n# Confidence model\nZ_MOD = 1.1     # Arbitrary modification to credibility coefficient","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"---\n# Tabular data preparation\n\n1. **Combination of observations** that occur in the same week for a given patient (keeps the average).\n\n2. **Smoothing of observations**: since measurements are imprecise (which is reflected in the fact that confidence values are floored at 70), we attempt here to eliminate some of the noise in the variation of measurements over neighboring weeks for a given patient. We hypothesize that each measurement is a normally-distributed random variable, and that series of measurements that are close to each other in time can be considered as realizations of multivariate normal distributions. This operation is done on *training* data, in an effort to help the model focus on signal rather than noise, but we keep validation data as is to protect our modelling of prediction error.\n\n3. **Interpolation between observations**: this step increases the amount of training data by linearly interpolating \"observations\" in weeks between given measurements.\n\n4. **Fold splitting**: since our end model is an ensemble, we randomly split the patients from the training set into groups that will be fixed throughout training and validation. While this creates an unenven number of observations for each fold, it ensures that the validation sets are composed exclusively of patients each component model has never seen, avoiding leakage."},{"metadata":{},"cell_type":"markdown","source":"## Combination of observations"},{"metadata":{"trusted":true},"cell_type":"code","source":"def combine_duplicates(df, FUN=np.mean):\n    \n    df = df.assign(patientWeeks = df.Patient + \"__\" + df.Weeks.astype(str))\n    table = df.patientWeeks.value_counts()\n    duplicates = table.loc[table > 1].index\n    subset = df.loc[df.patientWeeks.isin(duplicates)]\n    avgFVC = subset.groupby(['Patient', 'Weeks']).FVC.agg(FUN)\n    avgPct = subset.groupby(['Patient', 'Weeks']).Percent.agg(FUN)    \n    subset = subset.drop_duplicates(subset=['Patient', 'Weeks']\n                                   ).drop(labels=['FVC', 'Percent'], axis=1)\n    subset = subset.join(avgFVC, on=['Patient', 'Weeks']\n                        ).join(avgPct, on=['Patient', 'Weeks'])\n    df = pd.concat([df[~df.patientWeeks.isin(subset.patientWeeks)],\n                    subset]).sort_values(by=['Patient', 'Weeks'])\n    return df.drop(labels=['patientWeeks'], axis=1)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Smoothing of observations\n\nWe first compute a matrix of weekly relativities: if *n* is the total number of measurements and *d* is the maximum time lag/lead considered, creates a (*n*)x(2*d*+1) matrix, where each line represents, for a given patient at a given time *t0*, observations in the \\[*t0*-*d*, *t0*+*d*\\] range as a ratio of the observation at time *t0*. Since observations are not in fact available for every week, this step generates a NumPy Masked Array.\n\nFor each patient, we then consider each vector of relativities as a normally-distributed random vector **X** ~ N(***mu***, ***cov_matrix***), and seek the ***mu*** parameters that maximize the likelihood of our sample (we assume that covariance is fixed):\n1. We first compute the mean weekly relative values (as a baseline ***mu***) for the entire dataset.  \n\nThen, for each line of the matrix, we:\n2. Compute a covariance matrix (***cov_matrix***) for the weeks for which we have data;\n3. Search for a ***mu*** vector that maximizes the likelihood of the observations;\n4. Stack these optimized vectors in a matrix and return a weighted average of each column (giving more weight to rows that contain more data, to avoid giving undue influence to the edges of the time series) as the smoothed observations.\n\nThe calculation of the sample's likelihood contains a \"quick-and-dirty\" adjustment to account for the fact that some covariance matrices turn out to be negative definite. Several papers in the past two decades have proposed algorithms to find the closest positive definite approximation, but these do not seem to be implemented in Python functions yet, and wrapping the author's head around them would have required a disproportionate amount of effort (using absolute values in key places appear, graphically, to yield good-enough results for our current purpose)."},{"metadata":{"trusted":true},"cell_type":"code","source":"def weekly_relative_observations(data, delta):\n    \"\"\"\n    Outputs a matrix (NumPy Masked Array) of observed relativities.\n    \"\"\"\n    limit = delta * 2 + 1\n    positions = np.arange(-delta, delta + 1) # Relative weeks of interest\n    patients = data.Patient.unique()\n    obs_matrix = np.ma.masked_array([], mask=[])   # Observations matrix\n    for patient in patients:\n        subset = data.loc[data.Patient == patient]\n        weeks = np.array(subset.Weeks)\n        for i in range(len(weeks)):        \n            baseline = weeks[i]\n            rel_weeks = weeks - baseline\n            obs_exists = np.in1d(positions, rel_weeks) # Reverse mask\n            raw_obs_vec = positions.copy() * 0\n            relevant_weeks = weeks[np.in1d(rel_weeks, positions)]\n            if len(relevant_weeks) < 2:\n                continue\n            raw_obs_vec[obs_exists] = subset.loc[subset.Weeks.isin(relevant_weeks)].FVC\n            obs_vec = np.ma.masked_array([raw_obs_vec], mask=[~obs_exists])\n            if len(obs_matrix) == 0:\n                obs_matrix = obs_vec\n            else:\n                obs_matrix = np.ma.concatenate([obs_matrix, obs_vec], axis=0)  \n    return obs_matrix\n\n\ndef neg_log_likelihood_approx(mu, x, cov_matrix):\n    \"\"\"\n    Objective function to minimize, to smooth observations probabilistically.\n    \"\"\"\n    sq_mahalanobis = abs((x - mu).T.dot(np.linalg.inv(cov_matrix)).dot(x - mu))\n    return np.log(abs(np.linalg.det(cov_matrix))) + sq_mahalanobis\n\n\ndef probabilistic_smoothing(subset, obs_matrix, maxdelta=5):\n    \"\"\"\n    This function outputs a relativity vector, applicable to FVC and Percent \n    observations.\n    \"\"\"      \n    weeks = np.array(subset.Weeks)\n    FVC = np.array(subset.FVC)\n    target_weeks = np.arange(-maxdelta, maxdelta + 1)\n    mu = (obs_matrix / obs_matrix[:, maxdelta][:,np.newaxis]).mean(0)\n    n, m = len(weeks), max(weeks) - min(weeks) + 1\n    relativity_mat = np.zeros((n, m))     \n    for i in range(n):\n        week = weeks[i]\n        rel_weeks = weeks - week\n        relevant_obs = np.in1d(rel_weeks, target_weeks)\n        x = FVC[relevant_obs] / FVC[i]    \n        if len(x) <= 1:\n            relativity_mat[i, week - min(weeks)] = 1\n        else:\n            relevant_params = np.in1d(target_weeks, rel_weeks)\n            mu_mod = mu[relevant_params]\n            obs_matrix_mod = obs_matrix[:, relevant_params]\n            cov_mat = np.ma.cov(obs_matrix_mod, rowvar=False).data\n            mu_mod = minimize(neg_log_likelihood_approx, mu_mod,\n                              args=(x, cov_mat),\n                              method='Nelder-Mead')\n            relativity_mat[i, weeks[relevant_obs]-min(weeks)] = mu_mod.x\n    positive_rel = relativity_mat > 0\n    row_weights = positive_rel * positive_rel.sum(1)[:,np.newaxis]\n    keep_cols = relativity_mat.sum(0) > 0\n    return np.average(relativity_mat[:, keep_cols], axis=0,\n                      weights=row_weights[:, keep_cols])\n\n\ndef smooth_data(df, target_vars, maxrange=MAX_SMOOTHING_RANGE):\n    \"\"\"\n    Attempts to reduce the noise in each patient's time series data,\n    by smoothing variations in periods of frequent measurements. \n    \"\"\"\n    obs_matrix = weekly_relative_observations(df, maxrange)\n    for patient in df.Patient.unique():\n        target_loc = df.Patient == patient\n        rels = probabilistic_smoothing(df.loc[target_loc], obs_matrix, maxrange)\n        for target in target_vars:\n            df.at[target_loc, target] = df.loc[target_loc, target] / rels\n    return df","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Interpolation of observations"},{"metadata":{"trusted":true},"cell_type":"code","source":"def interpolate(df):\n    \"\"\"\n    Linearly interpolates the values of FVC and Percent between weeks with data\n    \"\"\"\n    ids = np.unique(df.Patient)\n    df_mod = pd.DataFrame()\n    for i in range(len(ids)):\n        subset = df.loc[df.Patient == ids[i]]\n        df_mod = pd.concat([df_mod, subset])\n        age, sex, smSt = subset.iloc[0].loc[['Age', 'Sex', 'SmokingStatus']]\n        for t in range(subset.shape[0]-1):\n            gap = subset.Weeks.iloc[t+1] - subset.Weeks.iloc[t]\n            if gap > 1:\n                base_Week, base_FVC, base_Pct = subset.iloc[t].loc[['Weeks', 'FVC', \n                                                                    'Percent']]\n                end_FVC, end_Pct = subset.iloc[t+1].loc[['FVC', 'Percent']]\n                for j in range(1, gap):\n                    new = pd.DataFrame({'Patient': ids[i],\n                                        'Weeks': base_Week + j,\n                                        'FVC': base_FVC + j / gap * (end_FVC - base_FVC),\n                                        'Percent' : base_Pct + j / gap * (end_Pct - base_Pct),\n                                        'Age': age, 'Sex': sex, 'SmokingStatus': smSt\n                                       }, index=[None])\n                    df_mod = pd.concat([df_mod, new], ignore_index=True)\n    return df_mod.sort_values(by=['Patient','Weeks']).reset_index(drop=True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Defining the tabular data container\n\n* Reads the training data CSV file;\n* Completes the aforementioned data preparation steps;\n* Splits the dataset into a specified number of folds with mutually exclusive validation sets;\n* Defines functions to pull training/validation data as needed."},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"scrolled":false},"cell_type":"code","source":"class CSVDataPrep():\n    \n    def __init__(self, data_path, nfolds):\n        data = pd.read_csv(data_path)\n        data = combine_duplicates(data)\n        self.raw_data = data\n        patients = np.array(data.Patient.unique())\n        self.plotter(self.raw_data, patients, True)\n        smooth = smooth_data(data.sort_values(by='Weeks'), ['FVC', 'Percent'])\n        self.plotter(smooth, patients, False)\n        self.data = []\n        for df in [data, smooth]:\n            # Saves unsmoothed and smoothed observations in positions 0 and 1 respectively\n            df = interpolate(df)\n            self.data.append(df.assign(targetFVC = df.FVC))\n        self.split_folds(data, nfolds)\n        plt.tight_layout()\n        plt.show()\n        \n        \n    def split_folds(self, data, nfolds):\n        nfolds = 2 if nfolds < 2 else int(nfolds)\n        patients = np.array(data.Patient.unique())\n        valid_size = int(np.round(len(patients) / nfolds))\n        self.train_sets, self.valid_sets = [], []\n        remaining = patients.copy()\n        for i in range(nfolds):\n            if i < nfolds - 1:\n                self.valid_sets.append(random.sample(list(remaining), valid_size))\n            else:\n                self.valid_sets.append(remaining)\n            self.train_sets.append(patients[~np.in1d(patients, self.valid_sets[i])])\n            remaining = remaining[~np.in1d(remaining, self.valid_sets[i])]\n        \n        \n    def pull(self, fold, raw=False):\n        \"\"\"\n        Pulls smoothed training and unsmoothed validation data for a single fold.\n        \n        If raw==True, omits interpolated weeks.\n        \"\"\"\n        fold = fold % len(self.train_sets)\n        data = self.raw_data\n        # Use smoothed data for training; unsmoothed data for validation\n        if not raw: data = self.data[1]\n        train = data.loc[data.Patient.isin(self.train_sets[fold])]\n        if not raw: data = self.data[0]\n        valid = data.loc[data.Patient.isin(self.valid_sets[fold])]\n        return train, valid\n    \n    \n    def plotter(self, dataset, patients, firstrun, per_row=4):\n        n = len(patients)\n        color = 'tab:cyan' if firstrun else 'tab:orange'\n        label = \"Raw\" if firstrun else \"Smooth\"\n        if firstrun:\n            fig, self.plots = plt.subplots(nrows=int(np.ceil(n/per_row)), ncols=per_row, \n                                           figsize=(12, n//per_row * 2.5))\n        for i in np.arange(n):\n            patient = patients[i]\n            subplot = self.plots[i//per_row, i%per_row]\n            subset = dataset.loc[dataset.Patient == patient]\n            x, y = np.array(subset.Weeks), np.array(subset.FVC)\n            subplot.plot(x, y, color=color, label=label)\n            subplot.set_xlabel(\"Weeks\")\n            subplot.set_ylabel(\"FVC\")\n            subplot.legend(loc=\"upper right\")\n        \n        \n        \ncsv_data = CSVDataPrep(DATA_PATH + \"train.csv\", N_FOLDS)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"---\n# Imaging data preparation\n\nHere we create a container to read, standardize and save imaging data. This class waits until a given patient's data is requested, then attempts to create a standardized approximation of the 3D scan data. It also creates a boolean array to let downstream parts of the algorithm know quickly if imaging data cannot be used for a given patient.\n\n* The dimensions of the 3D scan approximations are specified by the ***MAX_PLANES***, ***MAX_ROW*** and ***MAX_COL*** constants, defined at the top of the notebook.\n* Each 3D scan approximation takes a set of 2D scan approximations and create ***MAX_PLANES*** synthetic planes by aggregating the values of the ***N_NEIGHBOURS*** closest planes, using ***INVERSE_WEIGHT_FCT*** as the denominator when weighing the planes as a function of their distance to the synthetic plane.\n* Each 2D scan approximation samples down the .DCM pixel array using the appropriate sampling matrix.\n* The final result is a ((***MAX_PLANES***)\\*(***MAX_ROW***)\\*(***MAX_COL***))-dimension voxel vector. Admittedly, keeping a 3D structure and feeding it to a convolutional neural network might yield better results, but the author of this code is an AI novice who is still somewhat baffled by the process of selecting a good CNN architecture..."},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"class SamplingMatrix():\n    \"\"\"\n    Creates a matrix used for downsampling pixel arrays.\n    \"\"\"    \n    def __init__(self, source_size, target_size):\n        p = source_size / target_size\n        matrix = np.zeros((target_size, source_size))\n        for j in range(source_size):\n            for i in range(target_size):\n                if sum(matrix[:, j]) < 1:\n                    matrix[i, j] = min(1, p - sum(matrix[i,:]), 1 - sum(matrix[:,j]))\n                else:\n                    break\n        self.matrix = matrix / p\n\n        \n        \nclass SamplingMatrixContainer():\n    \"\"\"\n    Contains already-computed sampling matrices and computes a new one when required.\n    \"\"\"    \n    def __init__(self):\n        self.matrices = {}\n        \n        \n    def get_matrix(self, source_size, target_size):\n        if (source_size, target_size) not in self.matrices.keys():\n            self.matrices[(source_size, target_size)] = SamplingMatrix(source_size, \n                                                                       target_size)\n        return self.matrices[(source_size, target_size)].matrix\n\n\n    \nclass Standardized2DScan():\n    \"\"\"\n    This class loads and holds a pixel matrix from a given DICOM file. If the\n    matrix contains more rows or columns than the smallest image matrix in our\n    dataset, the __init__ method automatically shrinks it so that all matrices\n    used by our model have the same size.\n    \"\"\"    \n    def __init__(self, filepath, max_dim, sample_matrix_container):\n        with pydicom.dcmread(filepath) as f:\n            self.pixel_matrix = f.pixel_array\n            self.rows, self.cols = f.Rows, f.Columns\n            self.thickness, self.i = f.SliceThickness, f.InstanceNumber\n        del f\n        if self.rows > max_dim[0] or self.cols > max_dim[1]:\n            self.shrink_image(max_dim, sample_matrix_container)\n\n            \n    def shrink_image(self, max_dim, sample_matrix_container):\n        \"\"\"\n        Shrinks a pixel matrix to specified proportions.\n        \"\"\"\n        for i in range(2):\n            dim = self.pixel_matrix.shape[i]\n            if dim > max_dim[i]:\n                sampling_matrix = sample_matrix_container.get_matrix(dim, max_dim[i])\n                if i == 0: \n                    self.pixel_matrix = sampling_matrix.dot(self.pixel_matrix)\n                else: \n                    self.pixel_matrix = sampling_matrix.dot(self.pixel_matrix.T).T\n    \n    def get_vert_attr(self):\n        return self.thickness, self.i\n    \n    \n    \nclass Standardized3DScan():\n    \"\"\"\n    This class loads and holds the scan for a given patient, resizes each matrix if\n    necessary. It then creates a collection of equidistant synthetic planes along the\n    Z-axis by weighing the pixel values of the actual planes closest to each synthetic\n    plane, then stores the result in flat numpy array.\n    \"\"\"    \n    def __init__(self, patient, sample_matrix_container, trainpath=True,\n                 max_dim=(MAX_PLANES, MAX_ROW, MAX_COL), \n                 inv_weight_fct=INVERSE_WEIGHT_FCT, scaling_factor=PIXEL_VALUE_RANGE): \n        filepath = DATA_PATH + ('train' if trainpath else 'test') + '/' + patient + '/'\n        if os.path.isdir(filepath) == False:\n            self.voxel_vector = None\n            return None\n        files = np.array(os.listdir(filepath))\n        matrices = {}\n        vert_order = pd.DataFrame()\n        problem_files, self.error = [], False\n        for file in files:\n            try:\n                matrix = Standardized2DScan(filepath + file, max_dim[1:], \n                                            sample_matrix_container)\n                vert_order = vert_order.append({'file': file, 'location': matrix.i,\n                                                'thickness': matrix.thickness\n                                               }, ignore_index=True)\n                matrices[file] = matrix\n            except:\n                problem_files.append(file)\n                \n        if len(matrices) < max_dim[0]:\n            self.voxel_vector = None\n            return None\n        else:\n            # Deriving approximate .SliceLocation attributes\n            ordered = vert_order.sort_values(by='location')\n            ordered = ordered.assign(z = np.cumsum(ordered.thickness) - 0.5).sort_index()\n            z = np.array(ordered.z)\n            # Filter out problem files\n            files = files[~np.in1d(files, problem_files)]\n            ordered_array = np.concatenate([matrices[file].pixel_matrix \n                                            for file in files[z.argsort()]]\n                                          ).reshape(len(files), max_dim[1], max_dim[2])\n            z_matrix = self.z_weights(max_dim[0], z, ordered_array, inv_weight_fct)\n            voxel_array = np.zeros(np.prod(max_dim)).reshape(max_dim)\n            for i in range(max_dim[1]):\n                for j in range(max_dim[2]):\n                    voxel_array[:, i, j] = z_matrix.dot(ordered_array[:, i, j])\n            self.voxel_vector = voxel_array.flatten() / scaling_factor\n        del matrices\n        \n        \n    def z_weights(self, num_planes, z_vector, ordered_array, inv_weight_fct, \n                  n_closest_neighbours=N_NEIGHBOURS):\n        \"\"\"\n        Returns a matrix containing the weights of each actual plane for computing each \n        synthetic plane.\n        \"\"\"\n        if len(ordered_array.shape) != 3:\n            return False\n        z_range = max(z_vector) - min(z_vector)\n        n = ordered_array.shape[0]  # Original number of planes\n        weights_matrix = np.zeros((num_planes, n))\n        z_pos = np.array([(k+1)*(z_range) \n                          for k in range(num_planes)]) / (num_planes + 1) + min(z_vector)\n        for k in range(len(z_pos)):\n            z = z_pos[k]\n            distances = z_vector - z                                    \n            if np.any(distances == 0):\n                # No need to interpolate: we have an exact match\n                weights_matrix[k, np.where(distances == 0)[0][0]] = 1\n            else:\n                # Indices of closest points in distances vector\n                closest = np.absolute(distances).argsort()  \n                target_subset = closest[:n_closest_neighbours]\n                if np.all(z_vector[target_subset] > z) or np.all(z_vector[target_subset] < z):\n                    # We want to make sure we have at least one weighing point opposite \n                    # to the others\n                    # Smallest value greater than target z\n                    larger_index = min(np.where(np.sort(z_vector) > z)[0])  \n                    # Largest value lesser than target z\n                    smaller_index = max(np.where(np.sort(z_vector) < z)[0]) \n                    target_subset = np.concatenate([[larger_index], [smaller_index], \n                                                    target_subset[1:(n_closest_neighbours-1)]])\n                subset = distances[target_subset]  # Subset of distances of interest                \n                gross_weights = 1 / inv_weight_fct(np.absolute(subset))\n                weights_matrix[k, target_subset] = gross_weights / sum(gross_weights)\n        return weights_matrix\n\n    \n\nclass ScanDataContainer():\n    \"\"\"\n    Container for all standardized 3D arrays of scans. If a patient's scans have not yet been \n    processed, the get_voxel_vector() method automatically does it before returning them.\n    \"\"\"\n    def __init__(self, training=True, max_dim=(MAX_PLANES, MAX_ROW, MAX_COL), \n                 inv_weight_fct=lambda y: y):\n        self.filepath = DATA_PATH + ('train' if training else 'test') + '/'\n        self.max_dim = max_dim\n        self.inv_weight_fct = inv_weight_fct\n        self.images = {}\n        self.id_validation = {}\n        self.sample_matrix_container = SamplingMatrixContainer()\n        \n    def validate_availability(self, patients, trainpath=True):\n        for patient in patients:\n            if patient not in self.id_validation.keys():\n                validation = self.get_voxel_vector(patient, trainpath)\n                self.id_validation[patient] = validation is not None\n        return {p: self.id_validation[p] for p in patients}\n    \n    def get_voxel_vector(self, patient, trainpath=True):\n        if patient not in self.images.keys():\n            self.images[patient] = Standardized3DScan(patient, self.sample_matrix_container, \n                                                      trainpath=trainpath)\n        return self.images[patient].voxel_vector\n    \n    \n    \nimg_data =  ScanDataContainer()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"---\n# Prediction model\n\nWe create a pair of ensemble models, one that will use only tabular data, and one that will use both tabular and imaging data (both are contained in the **CustomEnsembleModel()** class, whose *predict()* method can call either model).\n\n## Image clustering\n\nThe size of the voxel vectors, given the computing power available and the small number of individuals in the training sample, calls for additional simplification of the imaging data. We attempt to drastically reduce dimensionality by running the K-means algorithm on our training set's imaging data, and transforming each observation's imaging data into a vector of euclidian distances to each centroid.\n\nFood for thought: perhaps a \"reverse-rendering\" approach might yield good dimensionality reduction while keeping valuable information about imaging features. This, however, lays rather significantly outside the author's current skillset."},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"class CustomClusterer():\n    \"\"\"\n    This class considers every training voxel vector instance as a centroid, and\n    transforms new image data as vectors of euclidian distances to each centroid.\n    \"\"\"\n    def __init__(self, img_data):\n        self.clusters = []\n        self.img_data = img_data\n        \n        \n    def fit(self, voxel_mat):\n        self.clusters = voxel_mat\n        \n        \n    def transform(self, voxel_mat):\n        result_matrix = np.array([])\n        for cluster in self.clusters:\n            distance = lambda vec1, vec2: np.linalg.norm(vec1 - vec2)\n            eucl_dist_vec = np.apply_along_axis(distance, 1, voxel_mat, cluster)\n            if result_matrix.shape[0] == 0:\n                result_matrix = eucl_dist_vec[:,np.newaxis]\n            else:\n                result_matrix = np.concatenate([result_matrix, \n                                                eucl_dist_vec[:,np.newaxis]\n                                               ], axis=1)\n        return result_matrix\n    \n    \n    \nclass ScanDataClusterer():\n    \"\"\"\n    This class repackages voxel data using the K-means algorithm OR the \n    CustomClusterer class. Outputs a DataFrame with the euclidian\n    distances between observations and each cluster.\n    \"\"\"    \n    def __init__(self, img_data, kernels=0):\n        self.img_data = img_data\n        self.n_kernel = kernels\n        if kernels <= 0:\n            self.kernelizer = CustomClusterer(img_data)\n        else:\n            self.kernelizer = KMeans(n_clusters=kernels)\n        self.scaler = StandardScaler()\n        \n    \n    def fit(self, patients):\n        voxel_mat = np.array([self.img_data.get_voxel_vector(p)\n                              for p in np.unique(patients)])\n        voxel_mat = self.scaler.fit_transform(voxel_mat)\n        self.kernelizer.fit(voxel_mat)\n        \n        \n    def fit_transform(self, data):\n        self.fit(np.array(data.Patient))\n        return self.transform(data)\n        \n        \n    def transform(self, newdata, trainpath=True):\n        patients = np.array(newdata.Patient)\n        voxel_mat = np.array([self.img_data.get_voxel_vector(p, trainpath=trainpath)\n                              for p in np.unique(patients)])\n        if len(voxel_mat.shape) == 1: \n            voxel_mat = voxel_mat.reshape(1, -1)    \n        voxel_mat = self.scaler.transform(voxel_mat)\n        voxel_df = pd.DataFrame(self.kernelizer.transform(voxel_mat))\n        \n        voxel_df.columns = [f'dist_k{k}' for k in range(len(voxel_df.columns))]\n        voxel_data = pd.DataFrame({'Patient': np.unique(patients)}).join(voxel_df)\n        indices = [voxel_data.loc[voxel_data.Patient == p].index[0] \n                   for p in patients]\n        return voxel_data.iloc[indices].reset_index(drop=True).drop('Patient', axis=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"markdown","source":"## Defining the predictive model\n\n* The **CustomEnsembleModel()** class contains 2 sets of component models for each fold of the training set: one that uses only tabular data, the other that also uses imaging data. Its predict method returns an average of the appropriate component models' predictions (bagging approach).\n* The **FibrosisModel()** class defines individual component models for the **CustomEnsembleModel()** class. It defines all the necessary preprocessing steps.\n* The **CustomModel()** class defines the neural networks that perform predictions for the **FibrosisModel()** class using data that has already been preprocessed."},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"class CustomModel():\n    \"\"\"\n    This class defines an underlying neural network for predictions, to be fed\n    fully preprocessed data.\n    \"\"\"\n    def __init__(self, shapes=[50,10]):\n        self.model = tf.keras.Sequential()\n        self.shapes = shapes\n        \n    def fit(self, X_train, y_train, eval_set, batch_size=50, epochs=2, \n            callbacks=[], verbose=0, dropout=0):\n        self.model.add(tf.keras.Input(shape=(X_train.shape[1],)))\n        for shape in self.shapes:\n            self.model.add(tf.keras.layers.Dropout(dropout, input_shape=(shape,)))\n            self.model.add(tf.keras.layers.Dense(shape, activation=\"elu\"))\n        self.model.add(tf.keras.layers.Dense(1))     \n        optimizer = tf.keras.optimizers.Adam()\n        self.model.compile(optimizer = optimizer,\n                           loss = tf.keras.losses.MeanSquaredLogarithmicError(),\n                           metrics = [tf.keras.metrics.MeanSquaredLogarithmicError()],)\n        history = self.model.fit(X_train, y_train, batch_size=batch_size, epochs=epochs,\n                                 validation_data=eval_set[0], #callbacks=callbacks, \n                                 verbose=verbose)\n        self.best_score = self.model.evaluate(*eval_set[0], batch_size=200)\n        \n    def predict(self, X):\n        return self.model.predict(X)\n\n        \n    def summary(self):\n        self.model.summary()\n\n        \n        \nclass FibrosisModel():\n    \"\"\"\n    This class encapsulates the underlying model as well as all preprocessing objects\n    and methods.\n    \"\"\"\n    def __init__(self, model_type, img_data=None, kernels=None, \n                 y=TARGET_VAR, **kwargs):\n        self.y = y\n        self.model = model_type(**kwargs)\n        self.img_data = img_data\n        if img_data is not None and kernels is not None:\n            self.kernelizer = ScanDataClusterer(img_data, kernels=kernels)\n        else:\n            self.kernelizer = None\n        self.cat_encoder = OneHotEncoder(handle_unknown='ignore', sparse=False)\n        self.scaler = StandardScaler()\n        self.firstrun = True\n        \n        \n    def sampler(self, data, sample_size):\n        data = data.reset_index(drop=True)\n        output = pd.DataFrame()\n        big_urn = np.array(data.index)\n        for i in range(sample_size):\n            choice = random.choice(big_urn)\n            draw = data.iloc[choice]\n            patient = draw.Patient\n            curWeek = draw.Weeks            \n            small_urn = np.array(data.loc[(data.Patient == patient) &\n                                          (data.Weeks != curWeek)].index)\n            choice = random.choice(small_urn)\n            target = data.iloc[choice]\n            draw.at[self.y] = target.loc[self.y]\n            output = output.append({**draw, 'targetWeek':target.Weeks\n                                   }, ignore_index=True)\n        return output.reset_index(drop=True)\n    \n    \n    def split_from_target(self, data):\n        return data.loc[:,~data.columns.isin(['Patient', self.y])], data.loc[:, self.y]\n    \n    \n    def preprocess(self, data, trainpath=True, img_data=None,\n                   cat_vars=['Sex','SmokingStatus'], scale_vars=['Percent','Age']):  \n        data = data.reset_index(drop=True)\n        if self.firstrun:\n            scale_cols = pd.DataFrame(self.scaler.fit_transform(data[scale_vars]),\n                                      index=data.index)\n            cat_cols = pd.DataFrame(self.cat_encoder.fit_transform(data[cat_vars]), \n                                    index=data.index)\n            if self.kernelizer is not None: \n                img_cols = self.kernelizer.fit_transform(data)\n            self.firstrun = False\n        else:\n            scale_cols = pd.DataFrame(self.scaler.transform(data[scale_vars]), \n                                      index=data.index)\n            cat_cols = pd.DataFrame(self.cat_encoder.transform(data[cat_vars]), \n                                    index=data.index) \n            if self.kernelizer is not None: \n                img_cols = self.kernelizer.transform(data, trainpath=trainpath)\n                \n        if self.img_data is not None and self.kernelizer is None:\n            img_cols = pd.DataFrame([self.img_data.get_voxel_vector(p, trainpath=trainpath)\n                                     / PIXEL_VALUE_RANGE\n                                     for p in np.array(data.Patient)])\n            img_cols.columns = [\"v\" + str(i) for i in range(img_cols.shape[1])]\n                \n        scale_cols.columns = scale_vars\n        cat_cols.columns = [s.replace(\" \",\"\").replace(\"-\",\"\") \n                            for s in np.concatenate(self.cat_encoder.categories_)]\n        data = data.drop(np.concatenate([cat_vars, scale_vars]), axis=1)\n        data = pd.concat([data, scale_cols, cat_cols], axis=1)\n        if self.img_data is not None:\n            data = data.join(img_cols)\n        return data.loc[:,np.sort(data.columns)]\n    \n    \n    def remove_invalid(self, data, trainpath=True):\n        if self.img_data is not None:\n            patients = np.array(data.Patient.unique())\n            availability = self.img_data.validate_availability(patients, trainpath)\n            patients = patients[[availability[p] for p in patients]]\n            return data.loc[data.Patient.isin(patients)]\n        else:\n            return data\n        \n        \n    def fit(self, train, valid, sample_size=1000, plot=False, early_stop=5, \n            verbose=False, callbacks=[], **kwargs):  \n        train, valid = self.remove_invalid(train), self.remove_invalid(valid)\n        train, valid = self.sampler(train, sample_size), self.sampler(valid, sample_size)\n        train, valid = self.preprocess(train), self.preprocess(valid)\n        X_train, y_train = self.split_from_target(train)\n        X_valid, y_valid = self.split_from_target(valid)\n        self.model.fit(X_train, y_train, eval_set=[(X_valid, y_valid)],\n                       verbose=verbose, #callbacks=callbacks, \n                       **kwargs)\n            \n    def predict(self, newdata, trainpath=True):\n        newdata = self.preprocess(newdata, trainpath)\n        newdata = newdata.drop(\"Patient\", axis=1)\n        if TARGET_VAR in newdata.columns:\n            newdata = newdata.drop(TARGET_VAR, axis=1)\n        return self.model.predict(newdata)\n    \n    \n    \nclass CustomEnsembleModel():\n    \"\"\"\n    This class creates an ensemble of models using k-folds validation\n    and averaging.\n    \"\"\"    \n    def __init__(self, csv_data, no_img_params, w_img_params, fit_params,\n                 img_data=img_data, nfolds=N_FOLDS):\n        self.csv_data = csv_data\n        self.img_data = img_data\n        self.nfolds = nfolds\n        self.models = {True: [], False: []}\n        self.model_params = {True: w_img_params,\n                             False: no_img_params}\n        self.fit_params = fit_params\n        self.score = {}\n        \n    \n    def evaluate(self, valid_set, with_imgs=False):\n        model = self.models[with_imgs]\n        rmse_vec = []\n        for i in range(self.nfolds):\n            X = model[i].remove_invalid(valid_set[i])\n            labels = np.array(X.loc[:,TARGET_VAR])\n            preds = model[i].predict(X)\n            if len(labels.shape) > 1:\n                labels = np.squeeze(labels)\n            if len(preds.shape) > 1:\n                preds = np.squeeze(preds)\n            rmse_vec.append(np.mean((labels - preds)**2)**0.5)\n        self.score[with_imgs] = np.mean(rmse_vec)\n        return self.score[with_imgs]\n        \n        \n    def fit(self, with_imgs=False, verbose=0, callbacks=[]):\n        model_set = self.models[with_imgs]\n        params = self.model_params[with_imgs]\n        img_data = self.img_data if with_imgs else None\n        \n        for i in range(self.nfolds):\n            print(\"\\nTraining model\", (i+1),\"/\",self.nfolds,\n                  (\"with\" if with_imgs else \"without\"),\"images\")\n            train, valid = self.csv_data.pull(i)\n            model = FibrosisModel(CustomModel, img_data=img_data, **params)\n            model.fit(train, valid, verbose=verbose, #callbacks=callbacks, \n                      **self.fit_params)\n            while model.model.best_score[0] > 0.1:\n                model = FibrosisModel(CustomModel, img_data=img_data, **params)\n                model.fit(train, valid, **self.fit_params)\n            model_set.append(model)\n            \n            \n    def predict(self, newdata, with_imgs=False, trainpath=True, folds=None):\n        pred_matrix = np.array([])\n        model_set = self.models[with_imgs]\n        if folds is None:\n            folds = range(len(model_set))\n        for fold in folds:\n            pred_vec = model_set[fold].predict(newdata, trainpath=trainpath)\n            if len(pred_vec.shape) == 1:\n                pred_vec = pred_vec[:, np.newaxis]\n            if len(pred_matrix.shape) == 1:\n                pred_matrix = pred_vec\n            else:\n                pred_matrix = np.concatenate([pred_matrix, pred_vec], axis=1)\n        try:\n            result = pred_matrix.mean(1)\n        except:\n            print(\"CustomEnsembleModel.predict() L171\")\n            print(pred_matrix.shape)\n            result = pred_matrix\n        return pred_matrix.mean(1)\n\n    \nmodel = CustomEnsembleModel(csv_data, no_img_params, w_img_params, fit_params)\nmodel.fit(False) # Train the version with only tabular data\nmodel.fit(True)  # Train the version that also uses imaging data","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"markdown","source":"---\n# Confidence model\n\nEstimates absolute prediction error by weighing the predictions of a set of linear models:\n* The different models predict absolute log errors as a function of a third-degree polynomial transformation of the time difference between the baseline week and the week for which a prediction is required.\n* We make an error model for the ensemble models and for each one of their component models, averaging the latter.\n* We hypothesize that the ensemble model (which is the one we ultimately use) will perform better on the test set than the component models do on their respective validation sets, but worse than the ensemble model performs on the training set.\n* The challenge is therefore to find a reasonable way of weighing the average of component errors against the ensemble error. We use an adaptation the Bühlmann-Straub credibility framework, which produces an interpolation coefficient ***z*** in the (0,1) range, used in actuarial science to assess the credibility of the experience of particular groups in a heterogenous portfolio. The value of ***z*** tends to 1 if we have many observations and/or a strong variance in the mean values of different groups; on the other hand, if there is little data to go by or if the experience tends to be stable across groups but highly variable within a given group, the value will tend to 0.\n* Let ***E*** and ***C*** be the error predicted, respectively, by the linear model trained from the known errors the ensemble model and the average of the linear models trained from the known errors of component models. We hypothesize that if there is strong heterogeneity between the component models' results, our overall predictions will be less reliable. Therefore, our final weighing will be:  \n**predicted error = *z* \\* *E* + (1-*z*) \\* *C***"},{"metadata":{"trusted":true},"cell_type":"code","source":"class ConfidenceModel():\n    \n    def __init__(self, prediction_models, tabular_data, img_data):\n        self.prediction_models = prediction_models\n        self.tabular_data = tabular_data\n        self.img_data = img_data\n        self.get_absolute_log_errors()\n        self.init_regressors()\n        self.compute_experience_params()\n        print(\"Done.\")\n        plt.show()\n    \n    \n    def get_absolute_log_errors(self, folds=N_FOLDS):\n        print(\"Computing known absolute log-errors...\")\n        data = self.tabular_data\n        # Prepare sets of predictions to make\n        pred_sets = []\n        whole_set = pd.DataFrame()\n        for fold in np.arange(folds):\n            _, valid = data.pull(fold, raw=True)\n            pred_set = pd.DataFrame()\n            for i, row in valid.iterrows():\n                subset = valid.loc[valid.Patient == row.Patient]\n                subset = subset.assign(targetWeek = subset.Weeks,\n                                       targetFVC = subset.FVC)\n                subset.Weeks = int(row.Weeks)\n                subset.FVC = float(row.FVC)\n                subset.Percent = float(row.Percent)\n                subset = subset.loc[subset.Weeks!=subset.targetWeek]\n                pred_set = pd.concat([pred_set, subset], ignore_index=True)\n            pred_sets.append(pred_set)     \n            whole_set = pd.concat([whole_set, pred_set], ignore_index=True)\n        pred_sets.append(whole_set)\n        # Compute absolute log errors\n        self.error_sets = {}\n        models = self.prediction_models\n        for img_context in [False, True]:\n            abs_log_error_set = []\n            for fold in np.arange(folds + 1):\n                pred_set = pred_sets[fold]\n                if img_context:\n                    with_imgs = self.img_data.validate_availability(np.array(pred_set.Patient))\n                    valid = [with_imgs[patient] for patient in pred_set.Patient]\n                    pred_set = pred_set.loc[valid]\n                preds = models.predict(pred_set, with_imgs=img_context,\n                                       folds=([fold] if fold < folds else None))\n                abs_log_errors = np.absolute(np.log(preds / pred_set.targetFVC))\n                ale_df = pd.DataFrame({'delta': np.array(pred_set.targetWeek - pred_set.Weeks),\n                                       'abs_log_error': abs_log_errors})\n                abs_log_error_set.append(ale_df)\n            self.error_sets[img_context] = abs_log_error_set\n            \n            \n    def init_regressors(self):\n        print(\"Initializing regressors...\")\n        self.regs = {}\n        self.limits = [0,1]\n        for img_context in [True, False]:\n            regs = []\n            for error_set in self.error_sets[img_context]:\n                reg = LinearRegression()\n                x = np.array(error_set.delta)\n                y = np.array(error_set.abs_log_error).reshape(-1,1)\n                if min(x) < self.limits[0]: self.limits[0] = min(x)\n                if max(x) >= self.limits[1]: self.limits[1] = max(x) + 1\n                x = np.stack([x, x**2, x**3], axis=-1)\n                reg.fit(x, y)\n                regs.append(reg)\n            self.regs[img_context] = regs\n            \n        \n    def compute_experience_params(self):\n        print(\"Computing Bulhmann-Straub credibility parameters...\")\n        limit = -self.limits[0]\n        deltas = np.delete(np.arange(*self.limits), limit)\n        self.z = {}\n        for img_context in [True,False]:\n            # Compute weights and average values by component/time-delta\n            weights = np.zeros((N_FOLDS, len(deltas)))\n            X = np.zeros((N_FOLDS, len(deltas)))\n            for fold in np.arange(N_FOLDS):\n                for delta in deltas:\n                    error_set = self.error_sets[img_context][fold]\n                    j = delta + limit + (-1 if delta>0 else 0)\n                    weights[fold, j] = sum(error_set.delta == delta)\n                    X[fold, j] = error_set[error_set.delta == delta].abs_log_error.mean()\n            X[np.isnan(X)] = 0\n            # Compute nonparametric Buhlmann-Straub quantities\n            Xww = np.average(X, weights=weights)\n            Xiw = np.average(X, weights=weights, axis=0)\n            Wiw = weights.sum(0)\n            s2 = (weights * (X - Xiw)**2).sum() / sum(Wiw-1)\n            a = sum(Wiw * (Xiw - Xww)**2) - (X.shape[1]-1) * s2\n            a = a * weights.sum() / (weights.sum()**2 - sum(Wiw**2))\n            z = Wiw / (Wiw + s2 / a)\n            z_df = pd.DataFrame({'delta': deltas, 'z': z})\n            self.z[img_context] = z_df\n            plt.plot(z_df.delta, z_df.z, color=('magenta' if img_context else 'cyan'),\n                     label=('Tabular + imaging data' if img_context else 'Tabular data'))\n        plt.xlabel(\"Weeks delta\")\n        plt.ylabel(\"Credibility coefficient\")\n        plt.legend(loc=\"center\")\n        \n        \n    def predict(self, start_week, img_avail=False, fwd=True, z_mod=Z_MOD):\n        regs, z = self.regs[img_avail], self.z[img_avail] # Get the right models     \n        # Time deltas of interest\n        mod = 1 if fwd else -1\n        weeks = np.arange(start_week, (MAX_TEST_WEEK if fwd else MIN_TEST_WEEK), mod) + mod\n        deltas = weeks - start_week  \n        # Credibility coefficients\n        coefs = np.array([(float(z.loc[z.delta==d].z) \n                           if d in np.array(z.delta) else 0)\n                          for d in deltas]) * z_mod\n        coefs = np.stack([coefs, 1 - coefs], axis=-1)\n        # Predictions\n        deltas = np.stack([deltas, deltas**2, deltas**3], axis=-1)\n        ens_pred = self.regs[img_avail][-1].predict(deltas)\n        comp_pred = np.concatenate([self.regs[img_avail][i].predict(deltas)\n                                    for i in np.arange(N_FOLDS)], axis=1).mean(1)\n        composite = np.concatenate([ens_pred, comp_pred.reshape(-1,1)], axis=1)\n        return np.exp((coefs * composite).sum(1)) - 1\n        \n            \nconfidence_model = ConfidenceModel(model, csv_data, img_data)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"---\n# Tying it all together...\n\nHere we pull the *test.csv* file and make predictions for all possible weeks, for each patient in the set, along with confidence values. Finally we save the result in as *submission.csv*."},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"def predict_all(model, data, img_data, confidence, trainpath=False,\n                lower=MIN_TEST_WEEK, upper=MAX_TEST_WEEK, default_conf=70):\n    \n    drop_at_the_end = [c for c in data.columns if c not in ['Weeks', 'Patient']]\n    weeks = np.arange(lower, upper+1)\n    output = pd.DataFrame()\n    data = data.assign(targetWeek = data.Weeks,\n                       targetFVC = data.FVC,\n                       Confidence = default_conf)\n    \n    # Visualization of results\n    per_row = 3\n    n = data.shape[0]\n    fig, axes = plt.subplots(nrows=int(np.ceil(n/per_row)), ncols=per_row,\n                             figsize=(12, n//per_row * 4))\n    \n    for idx, row in data.iterrows():\n        patient, start_week, start_FVC = row.Patient, row.Weeks, row.FVC\n        if len(output) > 0:\n            # To avoid any duplicates\n            if patient in np.array(output.Patient):\n                continue\n        \n        with_imgs = img_data.validate_availability([patient], trainpath=trainpath)\n        with_imgs = with_imgs[patient] # .validate_availability returns a dict\n            \n        past, future = weeks[weeks < start_week], weeks[weeks > start_week]\n        \n        pre_output = data.iloc[np.repeat(idx, len(weeks))].reset_index(drop=True)\n        pre_output.at[:, 'targetWeek'] = weeks\n        \n        for timerange in [past, future]:\n            if len(timerange) > 0:\n                if len(future) > 0:\n                    ascending = timerange[0] == future[0]\n                else:\n                    ascending = False\n                    \n                target_subset = pre_output.targetWeek.isin(timerange)\n                pred_subset = pre_output.loc[target_subset]\n                pred_subset = pred_subset.sort_values(by='targetWeek', ascending=ascending)\n                pred_subset = pred_subset.loc[:,~pre_output.columns.isin(['targetFVC', \n                                                                          'Confidence'])]\n                \n                # Take preds and confidence of imageless model in all cases\n                preds_no_img = model.predict(pred_subset, with_imgs=False, trainpath=trainpath)\n                conf_no_img = confidence.predict(start_week, fwd=ascending)\n                conf_no_img = np.maximum(1, np.absolute(conf_no_img * preds_no_img))\n                \n                # If imaging is available, average predictions from model with and without images,\n                # weighted by their respective confidence\n                if with_imgs:\n                    preds_w_img = model.predict(pred_subset, with_imgs=True, trainpath=trainpath)\n                    conf_w_img = confidence.predict(start_week, img_avail=True, fwd=ascending)\n                    conf_w_img = np.maximum(1, np.absolute(conf_w_img * preds_w_img))\n                    # Weighting\n                    weights = np.stack((conf_no_img, conf_w_img), axis=-1) ** -2\n                    weights = weights / weights.sum(1, keepdims=True)\n                    preds = (weights * np.stack((preds_no_img, preds_w_img), axis=-1)).sum(1)\n                    conf = (weights * np.stack((conf_no_img, conf_w_img), axis=-1)).sum(1) # Not quite rigorous but close enough\n                else:\n                    preds = preds_no_img\n                    conf = conf_no_img                    \n                        \n                if len(preds.shape) > 1:\n                    preds = preds[:,0]\n                if not ascending:\n                    preds = preds[::-1]\n                    conf = conf[::-1]\n                pre_output.at[target_subset,'targetFVC'] = preds\n                mod_conf = np.minimum(1000, conf) * 2 ** 0.5\n                pre_output.at[target_subset,'Confidence'] = mod_conf\n                \n                \n        \n        output = pd.concat([output, pre_output]).reset_index(drop=True)\n        \n        # Viz\n        if n <= per_row:\n            this_plot = axes[idx]\n        else:\n            this_plot = axes[idx//per_row, idx%per_row]\n        this_plot.set_xlim(lower, upper)\n        this_plot.set_ylim(0,4000)\n        this_plot.set_xlabel(\"Weeks\")\n        this_plot.set_ylabel(\"Predicted FVC\")\n        this_plot.set_title(patient)\n        x, y = pre_output.targetWeek, pre_output.targetFVC\n        yl, yu = y - pre_output.Confidence, y + pre_output.Confidence\n        this_plot.plot(x, yl, x, yu, color='cyan', linestyle='dashed')\n        this_plot.fill_between(x, yl, yu, where=yu>yl, facecolor='cyan', \n                               interpolate=True, alpha=0.5)\n        this_plot.plot(x, y, color='blue')\n        this_plot.axvline(x=start_week, color='k', linestyle='dashed')\n    \n    output = output.drop(drop_at_the_end, axis=1)\n    \n    plt.tight_layout()\n    plt.show()\n        \n    return output\n\n\ndef finalize_format(df):\n    df = df.assign(Patient_Week = df.Patient + '_' + df.targetWeek.astype(str))\n    df = df.rename(columns={'targetFVC':'FVC'})\n    return df.drop(['Patient', 'targetWeek', 'Weeks'], axis=1)\n\n\ntest_data = pd.read_csv(DATA_PATH + \"test.csv\")\noutput = predict_all(model, test_data, img_data, confidence_model)\nfinal = finalize_format(output)\nfinal.to_csv(\"submission.csv\", index=False)\nfinal.head()","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}